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1.0  Summary  of  the  Report 


The  two-year  study  involved  a  series  of  investigations  which  are 
summarized  below. 

1.1  Radar  Reflectivity  Studies 

The  factors  contributing  to  the  variability  of  radar  reflectivity  - 
water  content  relationships  are  investigated  using  a  combination  of 
numerical  simulations  using  the  3-D  cloud/mesoscale  model  and  field 
measurements.  It  was  fonnd  that  model  predictions  of  reflectivities  due 
to  ice  crystals  are  inaccurate  due  to  a  lack  of  information  on 
predicting  the  concentrations  of  ice  crystals  in  cloud.  In  contrast,  it 
was  found  that  model-predicted  PPI  and  RHI  reflectivities  compared 
favorably  to  measured  reflectivities  for  Florida  cumulus.  Raindrops  and 
graupel  particles  caused  the  measured  reflectivities;  ice  crystals 
contributed  negligibly.  The  model-predicted  reflectivities  were  deemed 
sufficiently  accurate,  from  the  Florida  comparison,  to  investigate  the 
causes  for  the  variability  in  water  content-radar  reflectivity 
rela  tionships . 

Aircraft  simulated  'penetrations'  of  cumulus  clouds  using  the  3-D 
model  revealed  regions  near  cloud  edge  where  water  contents  were  low  and 
reflectivities  were  high  due  to  the  presence  of  a  few  large  graupel 
particles.  This  finding  was  supported  by  an  analysis  of  aircraft 
penetrations  of  northeast  Colorado  hailstorms. 

Additional  analysis  of  these  hailstorm  data  revealed  a  relationship 

_3 

between  water-content  (M)  and  reflectivity  (Z)  of  M  (gm  )  =  0.065 
JjO.196  (m^m  ^).  The  exponent  of  0.196  varies  from  more  common  0.5  to 
1.0  values  due  to  the  presence  of  large  hail.  The  Marshall-Palmer  size 
distribution  used  in  the  3-D  model  causes  M  to  be  proportional  to  Z*’®. 
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Hence,  the  3-D  model  is  not  able  to  simulate  reflectivities  in 
hailstorms . 

1 .2  Aircraft  Icing  Studies 

The  3-D  cloud/mesoscale  model  was  adapted  to  simulate,  in  an 
orographic  cloud,  the  production  of  liquid  water  and  the  removal  by 
various  ice  crystal  processes.  To  adapt  the  model,  an  ice  crystal 
nucleation  parameterization,  an  ice  crystal  aggregation  parameterization 
and  an  ice  crystal  removal  parameterization  were  developed.  Preliminary 
simulations  have  revealed  an  interesting  relationship  between  ice 
crystal  nucleation  and  the  existence  of  supercooled  water.  The 
production  of  aggregates  also  affected  the  existence  of  liquid  water 
because  they  removed  large  numbers  of  ice  crystals  and  thus  affected  the 
overall  balance. 

1 .3  One-Dimensional  Cloud/Turbulence  Model  Investigations 

The  theory  of  the  one  dimensional  cloud/ turbulence  model,  which 
satisfactorily  simulates  many  features  of  the  cloud-capped  atmospheric 
boundary  layer,  was  extended  to  include  ice  processes.  As  a  result  of 
this  investigation,  a  simplified  model  was  being  developed. 

1 .4  Hydrostatic  Mesoscale  Model  Investigations 

The  non-hydrostatic  3D  cloud/mesoscale  model  was  successfully 
adapted  to  a  hydrostatic  mesoscale  model.  The  hydrostatic  mesoscale 
model  provides  greater  computational  efficiency  for  large-scale 
simulations.  A  data  assimilation  and  analysis  package  was  developed  to 
initialize  the  hydrostatic  model  with  standard  meteorological 


measurements . 
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1 • 5  Investigations  of  Wind  Shears  Prodnced  bv  Precipitating 
Convective  Clonda 

Strong,  transient,  low-level  wind  shears  were  analyzed  from  data 
collected  in  South  Park,  Colorado  and  Niles  City,  Montana  thunderstorms. 
Model  simulations  based  on  the  data,  reveal  similar  wind  shear  and  the 
causal  processes.  The  shears  are  primarily  associated  with  thunderstorm 
outflow  regions  which  were  driven  by  precipitation. 

1.6  Recent  Change  in  the  3-D  Cloud/Mesoscale  Model  Sent  to 
Kirtland  AFB. 

Recent  changes,  which  improved  the  current  model  formulation  are 


detailed . 
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2.0  Introduction 

The  estimation  of  water  contents  of  clouds  using  conventional  radar 
has  become  important,  recently,  because  zones  of  high  liquid  water 
contents  may  result  in  the  destruction  of  aerospace  vehicles  upon  re¬ 
entry  into  the  troposphere.  Empirical  relationships  between  the  liquid 
(or  ice)  water  content  of  precipitation  (M)  and  radar  reflectivity  (Z) 
generally  have  been  derived  for  clond  systems  with  one  predominant 
precipitation  phase  (eg.,  either  ice  or  water).  Recently,  Plank  et  al . 
(1980)  have  presented  empirical  M-Z  relationships  for  ice  hydrometeors 
using  data  collected  from  new  electro-optic  sensors  (Knollenberg,  1970). 

In  this  study  we  investigated  M-Z  relationships  in  mixed-phase 
clouds.  We  employed  both  measured  M-Z  relationships  and  relationships 
predicted  using  the  three-dimensional  cloud/mesoscale  numerical 
simulation  of  Tripoli  and  Cotton  (1982).  A  unique  feature  of  this 
research  was  the  coupling  of  measured  and  simulated  M-Z  relationships. 
Using  this  approach,  it  was  found  that  much  of  the  variability  in  the 
relationships  can  be  explained. 

2.1  Obiectives 

The  initial  objectives  of  the  proposal  (for  the  period  November 
1981  through  March  1983)  were  as  follows: 

•  To  compare  model-predicted  liquid  water  (M)  -  radar  reflectivity 
(Z)  relationships  with  measured  M-Z  relationships  from  Florida, 
Colorado  and  Montana  cumulus  clouds. 

*  To  determine  the  dominant  factors  contributing  to  the 
variability  in  M-Z  relationships  using  model-predicted  and 
measured  data. 

The  final  objectives  of  the  project  (for  the  period  November  1982  to 
September  1983)  were  as  follows: 
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*  Use  the  3-D  cloud/mesoscale  model  to  define  potential  regions 
for  aircraft  icing. 

*  Determine  the  feasibility  of  adapting  the  existing  1-D 
cloud/ turbulence  model  to  include  ice  processes  and  modify  the 
model  accordingly. 

*  Expand  the  CSD  3-D  cloud  model  into  a  hydrostatic  mesoscale 
model  to  study  regional  cloud  formations  with  application  to 
predictions  of  EM  energy  propagation,  low-level  wind  fields 
atmospheric  transport,  etc. 

*  Investigate  causes  of  extreme  wind  shears  in  the  atmosphere 

*  Descriptions  of  the  updates  to  the  CSD  cloud/mesoscale  mode 
which  is  currently  installed  on  the  Kirtland  AFB  computer  fo*. 


AFGL/LYC  use. 


6 


3.0  Radar  Reflectivity  Studies 

3 .1  Model  Predicted  Radar  Reflectivity 

As  our  first  step  in  analyzing  the  variability  of  liquid  water  (M) 
versus  radar  reflectivity  (Z) ,  we  have  introduced  algorithms  in  the  3D 
model  to  analyze  and  display  radar  reflectivity.  The  details  of  the 
algorithms  are  given  in  the  Appendix  (section  3.6).  The  standard 
definition  of  radar  reflectivity  is  of  the  form 


Z 


I 

vol 


(3.1) 


where  n^  represents  the  concentration  of  scattering  elements  per  unit 
volume  of  diameter  D^.  This  definition,  however,  does  not  include  any 
effects  of  the  phase  of  the  scattering  elements  which  modulate  the 
amplitude  of  the  complex  index  of  refraction. 

Therefore,  we  defined  a  complex  index  of  refraction-weighted 
reflectivity  Z  such  that 


Z  = 


III2  X 

vol 


(3.2) 


where  we  assume  |kp  is  0.93  x  10*^  for  liquid  particles  and  0.19  x  10*^ 
for  ice  particles.  Ice  particles  in  the  melting  zone  are  assumed  to  be 

i  i2  12  12 

water-coated  and  therefore  have  a  |k{  of  0.93  x  10  .  (The  10  is  a 

-3  ~ 

unit  conversion  factor  such  that  if  D.  is  in  cm  and  n.  in  cm  ,  Z  is  in 

l  l 


By  using  Z  we  can  examine  the  effects  of  particle  size 
distributions  and  the  phase  of  precipitation  elements  on  the  variability 


of  M  vs.  Z  relationships. 

In  the  3D  model  analysis  we  introduced  algorithms  for  computing  the 
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modified  reflectivity  for  raindrops  Z  ,  graopel  particles  Z  ,  ice 

r  S 

crystals  Z^,  and  aggregates  Z^.  Under  the  assumed  particle  size 
distribution  functions  in  the  model,  the  reflectivities  for  the 
individnal  components  are  as  follows. 

For  raindrops. 


Z 

r 


2 6  r<7)  r3 

_ m 

8n  p 

w 


r 

r 


1.7xl015  R  p  r 
m  o  r 


(3.3) 


where  R  is  a  characteristic  drop  radius  for  the  distribution,  p  is  the 

n  w 

density  of  water,  p  index  of  refraction  and  r  is  the  mixing  ratio  of 
o  r 

raindrops . 

For  graupel. 


Z 

g 


rmp 

- e-ai  |kJii(l+H(T-273.16)*3.9)r 


np 


g 


g 


11  ,  p  r 

4. 35x10^(1+^1-273. i6)*3.9)D*  -e-* 


mg  p 


g 


(3.4) 


where  p  is  the  density  of  a  graupel  particle,  D  is  the  characteristic 

g  <ng 

diameter  of  graupel,  |k^|^  is  the  ice-phase  complex  index  of  refraction, 

T  is  temperature,  H(x)  is  the  heaviside-step  function  and  r  is  the 

g 

mixing  ratio  of  graupel  particles.  Note,  the  formulation  differs  from 
the  first  quarterly  report  due  to  recent  reformulation  of  the  assumed 
graupel  distribution  from  the  assumption  of  constant  concentration  to 
the  assumption  of  constant  mean  diameter. 


For  ice  crystals 
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*.  ■  i‘J2  \  A 


(3.5) 


where  N.  is  the  ice  crystal  concentration  and  d.  is  the  diameter  of  the 
1  l 


ice  crystals. 


For  aggregates. 


r<7)D3 

Z  =  - —  p  Ik.l  [l+3.9H(T-273.16)]r 

a  np  o'  t'  a 

«g 


=  2.9xl015[l+3.9H(T-273.l6)]D3*6p  r 

ma  o  a 


(3.6) 


where  D  is  the  characteristic  diameter  of  aggregates,  |k|^  is  the 
na 

ice-phase  complex  index  of  refraction,  r&  is  the  mixing  ratio  of 

*K8regate  particles  and  p  is  the  aggregate  density  given  by 

*8 


p  =  0.015  D°‘6 

ag  ma 


(3.7) 


The  total  reflectivity  factor  Z^.  was  assumed  to  be 


Z  =  Z  +  Z.  +  Z 
T  g  i  a 


(3.8) 


This  can  be  compared  with  the  total  mass-density  of  all  precipitation 
elements  (M^.)  given  by: 


M_=p  (r  +  r  +  r ,  +  r  ) 
Tor  g  i  a 


(3.9) 


3 . 2  Simulations  of  a  South  Park,  Colorado  Thunderstorm 
The  analysis  procedure  was  applied  to  data  obtained  from  the 
simulation  of  a  quasi-steady  storm  observed  during  the  1977  South  Park 


Area  Canales  Experiment  (SPACE).  The  3D  model  was  used  to  predict 


rf,  r^  and  r^  and  the  corresponding  reflectivities.  Figure  3.1a,b.c 

show  the  fields  of  reflectivity  for  rainwater,  graupel  and  ice  crystals, 

respectively  in  units  of  dBZ,  i.e.,  10  l°gj0  (Z).  Figure  3.2  shows  the 

resultant  total  field  of  reflectivity  or  Z  in  dBZ.  The  total 

P 

reflectivity  field  shows  the  combined  influence  of  graupel  particles  in 
the  main  updraft  region  of  the  storm  and  ice  crystals  in  the  outflow 
region  at  about  the  6.5  km  level  in  the  rear  quadrant  of  the  simulated 
storm . 

These  illustrations  (Fig.  3.1  and  3.2)  demonstrate  some  of  the 
problems  encountered  when  trying  to  use  the  3D  model  to  predict  radar 
reflectivity.  Clearly  the  peak  55  dBZ  obtained  for  graupel  in  Fig.  3.1b 
is  realistic  and,  in  fact,  compares  well  with  radar  measurements  made  on 
that  day.  The  rain  radar  reflectivities  in  Fig.  3.1a  also  look 
realistic  although  their  peaks  are  confined  to  very  close  to  the  ground. 
The  ice  crystal  reflectivities,  however,  demonstrate  a  major  problem 
with  the  reflectivity  predictions.  As  we  can  see  in  Fig.  3.1c,  the 
reflectivity  generally  increases  downward  until  a  sharp  cut  off  at  the 
melting  zone  where  we  assume  the  crystals  melt  simultaneously.  This 
reduction  in  reflectivity  is  becanse  we  are  diagnosing  the  crystal 
concentration  from  the  Fletcher  temperature-concentration  relationship 
which  gives  approximately  a  one  order  of  magnitude  decrease  in 
concentration  for  each  4°C  rise  in  temperature.  Therefore,  ice  crystal 
populations  produced  near  -20°C  which  settle  to  the  melting  zone 
decrease  strongly  in  number  and  therefore  their  diagnosed  size  increases 
strongly.  Since  Z^  is  proportional  to  D^,  a  very  strong  increase  in  dBZ 
is  predicted.  So  strong  is  this  effect,  that  the  Z^  actually  completely 
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Figure  3.1b:  Same  as  3.1a  except  for  ice  crystals. 
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Figure  3.1c:  Same  as  3.1a  eacept  for  graupel. 


doainates  the  total  reflectivities  just  above  the  melting  zone. 


We  thought  this  problem  would  be  overcome  with  the  prediction  of 
ice  crystal  concentration,  since  large  numbers  of  crystals  higher  in  the 
cloud  would  lower  the  mean  sizes  below.  The  implementation  of  a  crystal 
prediction  has  been  accomplished  in  our  orographic  cloud  study.  The 
results  of  ice  crystals  reflectivities  produced  in  that  study  are  shown 
in  Figure  3.2a.  In  this  wintertime  case,  the  freezing  level  is  near  the 
gronnd.  Although  the  dominating  tendency  for  Z.  to  be  simply  a  function 
of  temperature  is  absent,  we  now  find  the  existence  of  random  pockets  of 
unrealistically  high  Z values  coexisting  with  realistic  Z.  values.  The 
random  pockets  occur  most  often  on  the  cloud  boundaries.  We  can  explain 
these  results  as  follows:  Unlike  all  other  quantities  predicted  in  the 
model,  is  poorly  behaved.  Whereas  r^  may  vary  by  one  order  of 
magnitude  in  the  cloudy  region,  may  vary  by  5-10  orders  of  magnitude. 
Consequently,  the  meaningful  prediction  of  transport  and  diffusion  of  N 
on  a  finite  difference  grid  will  be  much  less  accurate  than  for  the 

prediction  of  r.  especially  in  regions  where  N  is  changing  sharply. 

1  i 

Such  regions  are  most  often  at  the  cloud  boundaries.  It,  therefore, 
appears  that  the  prediction  of  improves  Z^  predictions  in  the  cloud 
interior,  but  leads  to  unacceptably  noisy  and  large  dBZ  patterns  near 
the  boundary  of  the  ice  crystal  region.  These  problems  with  ice  crystal 
concentration  predictions  are  less  severe  for  the  prediction  of  crystal 
mixing  ratios  since  the  noisy  regions  usually  have  very  small  ice 
contents.  However,  the  predictions  will  have  left  us  with  an  unreliable 
M-Z  relationship  for  ice  crystals,  near  the  cloud  boundaries,  although 
realistic  predictions  appear  possible  in  the  cloud  interior. 
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3.3  Measured  and  Predicted  Radar  Ref lectivities  for  Florida 
Cumulus 

The  measured  radar  reflectivities  were  obtained  from  the  University 
of  Miami  Doppler  radar  courtesy  of  John  Cunning  (1981,  personal 
communication) .  Cunning  reported  the  radar  was  simultaneously  operated 
with  the  similar  NCAR  CP-4  radar  (prime  calibration,  signal  generation) 
and  the  UM  radar  values  were  adjusted  to  the  CP-4  values.  Cunning  et 
al .  (1979)  presents  a  brief  discussion  of  the  visual  appearance  of  the 
investigated  clouds  on  23  August  1975  as  well  as  the  PPI  radar 
characteristics.  Here,  we  present  the  PPI  reflectivity  values  with 
superimposed  flight  tracks  of  the  sampling  aircraft  (NOAA-C130)  in  Fig. 
3.3. 

The  measured  reflectivity  (Z)  values  were  obtained  from  Fig.  3.3 
using  the  following  procedure.  Particle  spectra  shown  in  Fig.  3.4  were 
derived  from  the  entire  aluminum  foil  sample  for  a  penetration  of  the 
C-130  aircraft,  hence  producing  an  average  spectrum  and  water  content 
(M)  for  a  particular  penetration.  Consequently,  it  was  not  possible  to 
correlate  "M"  values  with  "Z"  values  as  they  varied  along  a 
penetration  path.  Instead,  area-averaged  ”Z"  values  along  a 
penetration  path  were  determined  from  Fig.  3.3.  The  results  are  listed 
in  Table  3.1.  It  can  be  seen  by  inspecting  Fig.  3.3  for  peak  "Z" 
values  and  Table  3.1  for  area-averaged  "Z"  values,  that  the  peak 
values  ranged  from  30  to  50  dBZ  while  the  average  values  range  from  18 
to  30  dBZ. 

The  calculated  and  predicted  radar  reflectivity  values  were 
obtained  from,  respectively,  the  particle  spectra  in  Fig.  3.4  and  from 
the  three-dimensional  numerical  simulations  of  Levy  (1982). 
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Figure  3.3:  Radar  reflectivity  values  (dBZ)  measured  on  25  August  1975 
in  South  Florida.  The  initial  contour  corresponds  to  15 
dBZ,  with  5  dBZ  increments  thereafter.  The  aircraft 
penetrations  are  given  by  the  dashed  lines  with 
corresponding  beginning  and  ending  times.  The  PPI  radar 
images  are  from  the  6  km  level  and  the  approximate  scan 
times  are  given. 
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Figure  3.4e:  Size  distributions  of  ice  and  liquid  particles  at  the 
-10°C  level  in  a  Florida  cumulus  tower,  25  August  1975, 
143443  to  143507  EDT.  The  undetermined  particles  are 
either  ice  or  water  or  both.  The  water  contents  derived 
from  the  size  distributions  are  given.  The  size 
distributiona  and  water  contents  are  accurate  to  +  50%, 
based  on  the  data  redaction  strategy  discussed  in  the 
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Figure  3.4b:  Particle  sire  distributions  and  water  contents  from 
adjacent  tower,  25  August  1975,  143752-143805  EDT. 
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The  reflectivities  were  calculated  from  the  particle  size  spectra 


using  the  following  expression: 


n 

z  =  r 

i-l 


N.D6 
1  1 


,  6  3 . 

[  mm  m  J 


(3.10) 


where  is  the  number  of  particles  in  diameter  interval  i  to  i  +  1  and 
is  the  melted  diameter  of  the  particles  in  the  interval  defined  by 


the  geometric  diameter  D  =  ,  D.  *  D 

g  M  1 


i+l 


Following  the  procedure 


developed  in  Section  3.1,  Eq .  (3.10)  becomes 


! k | 2  I  N.D6 
i=l  1  1 


(3.11) 


where  Z  is  an  index-of-refraction  weighted  reflectivity.  he  assume  Ul‘ 
is  0.93  for  liquid  particles  and  0.19  for  ice  particles.  Finally,  the 
total  reflectivity  value  Z^  is  defined  as 


Z  =  Z  4  Z  4  Z 

t  r  g  i 


(3.12) 


The  reflectivities  (Z  ,  Z  ,  Z.,  and  Z  )  calculated  from  the  particle 

T  g  1  t 

spectra  are  given  in  Table  3.1  in  terms  of  dBZ  (10  log  Z  -  dBZ) . 

Three  features  are  apparent  in  Table  3.1  between  the  radar-measured 
reflectivities  and  the  reflectivities  calculated  from  the  particle 
spectra.  First,  the  measured  and  calculated  total  reflectivities  track: 
25  -  30  vs.  14.6  dBZ  (Cell  A),  18  -  23  vs.  7.3  dBZ  (Cell  B)  and  25  vs. 
20.5  dBZ  (Cell  C) .  Second,  the  measured  values  are  consistently  greater 
than  the  calculated  values.  Third,  the  graupel  water  usually 
contributes  the  greatest  to  the  Z^  values. 

The  model  reflectivities  were  obtained  from  the  three-dimensional 
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numerics’  simulations  of  Levy  (1982)  following  the  area- averaging 

procedure  used  with  the  measured  PPI  reflectivities.  The  results  of  the 

three-dimensional  simulations  were  in  the  form  of  PPI  reflectivity  plots 

at  the  -10°C  level  (6.3  km,  aircraft  penetration  level)  as  a  function  of 

time  after  initiation  of  convection.  The  Z  ,  Z  ,  Z.  and  Z  values 

r  g  l  t 

derived  from  the  plots  are  listed  in  Table  3.1  in  terms  of  dBZ. 

Two  features  are  apparent  in  Table  3.1  between  the  radar-measured 

reflectivities  and  the  reflectivities  calculated  with  the  model.  First, 

the  measured  values  range  between  IS  and  30  dBZ  and  the  calculated 

values  range  between  17  and  23  dBZ.  This  result  demonstrates  a 

remarkable  agreement  between  the  measurements  and  predictions.  Second. 

the  model-predicted  Z  values  contribute  more  to  Z,  than  the  Z  values. 

r  t  g 

In  contrast,  the  Z  values  calculated  from  the  particle  spectra 

8 

contributed  more  to  Z  than  the  Z  values.  It  is  not  possible  from 

t  r  * 

these  conflicting  results  to  determine  which  type  of  precipitation 

contributes  the  most  to  the  reflectivity  values.  Nevertheless,  it  is 

clear  that  Z  and  Z  values  are  significant  contributors.  Additional 
r  g 

data  sets  are  required  to  determine  the  major  contribution.  This  was 
accomplished  in  Section  3.4  using  NHRE  data. 

Vertical  cross-sections  of  the  measured  radar  reflectivity  were 
also  compared  to  model-predicted  radar  cross-sections.  Measured 
reflectivities  shown  in  Fig.  3.5a-d  are  contoured  in  intervals  of  10  dBZ 
up  to  40  dBZ  and  at  5  dBZ  intervals  thereafter.  The  model-predicted 
reflectivities  in  Fig.  3.5e  and  3.5f  are  contoured  at  uniform  10  dBZ 
intervals.  In  general  the  model-predicted  reflectivities  are  the  same 
as  observed  and  the  general  cell  evolution  is  similar.  The  upper  level 
maximum  of  Z  predicted  by  the  model  corresponds  to  the  early  stages  of 

y 


graupel  formation.  Such  a  maximum  was  not  measured,  however. 


In  summary,  the  measured  PPI  reflectivity  values  were  consistently 

greater  than  the  reflectivities  calculated  from  the  particle  spectra. 

In  contrast,  the  measured  reflectivities  (PPI  and  RH 1 )  were  remarkably 

consistent  with  the  reflectivities  resulting  from  the  three-dimensional 

simulation.  The  Z  and  Z  values  contributed  to  the  Z  values.  It  was 
r  g  t 

not  possible  to  determine  the  dominant  precipitation  type  because  the  Z 

and  Z  values  resulting  from  the  particle  measurements  and  three- 
g 

dimensional  simulations  were  not  consistent. 

There  are  two  possible  explanations  for  these  discrepancies  between 
the  measured  and  calculated  reflectivities.  First,  the  radar-derived 
reflectivities  may  be  overestimated  because  the  reflectivities  are 
presumably  caused  by  liquid  particles.  We  have  shown  that  graupel  par¬ 
ticles  contribute  significantly  to  the  reflectivities.  Second,  the 
measure  of  rain  water  and  graupel  water  contents  may  have  been  underes¬ 
timated  due  to  limited  sample  volumes  leading  to  the  underestimated  Z 
values  from  the  particle  spectra.  Nevertheless,  differences  between  the 
calculations,  predictions  and  the  measurement  were  on  the  order  of  the 
uncertainties  of  the  measurements.  Consequently,  the  model- predicted 
reflectivities  are  suitable  for  investigating  the  causes  for  the  varia¬ 
bility  in  water  content-radar  reflectivity  relationships. 

3 .4  Investigations  of  Variations  in  Water  Content  (M)  -  Radar 
Reflectivity  (Z)  Relationships 
3.4.1  Model  investigations 

The  model  investigation  of  M-Z  variations  first  focused  on  investi¬ 
gating  ’simulated*  aircraft  penetrations.  The  results  from  the  3-D 
model  calculations  of  the  Florida  cumulus  were  used  to  construct  three 


'simulated'  aircraft  penetrations  at  the  6.3  km  (-10  C)  level  (see  Fig. 
3.6).  It  can  be  seen  in  Figure  3.6  that  the  simulations  consisted  of 
rain  and  graupel  water  contents  (p^  and  p^) ,  radar  reflectivities 
(lOlogZ'),  and  the  ratio  of  lOlogZ*  and  (p^  +  p^)  as  functions  of  dis¬ 
tance  in  the  cloud  [Z'  is  the  index  of  refraction-weighted  reflectivity, 
equal  to  Z  x  0.93  for  liquid  and  Z  X  0.19  for  ice].  The  latter  ratio 
was  calculated  to  evaluate  any  region  where  the  total  reflectivity 
behaved  anomalously  with  respect  to  the  total  mass.  The  significant 
result  from  these  simulations  are  the  regions  of  high  reflectivities  and 
low  water  contents  near  the  edge  of  the  cloud  caused  by  a  relatively  few 
large  particles.  The  water  contents  in  these  regions  are  less  than  10% 
of  the  maximum  water  contents  in  the  cloud.  Consequently,  the  anomalous 
regions  contribute  little  to  the  total  water  content  of  the  cloud.  But, 
at  low  water  contents,  a  sufficiently  large  number  of  large  particles 

are  produced  through  the  p^  <*  N(D)  parameterization  to  affect  the 
6  3 

reflectivity  («  D  )  but  not  the  water  content  («  D  ) .  Further,  in  the 
high  reflectivity  regions  in  the  center  of  the  cloud,  the  ratios  are 
constant  with  distance  in  the  cloud.  These  results  were  investigated 
further  in  the  next  section  to  determine  their  realism  based  on  compari¬ 
son  with  measurements  in  cumulus  clouds.  It  wll  be  shown  that  the  model 
predicted  "anamalous  M  vs.  2"  regions  are  physically  real. 

3.4.2  NHRE  case  study:  M  «  Z°‘2 
The  case  selected  for  investigation  was  a  hailstorm  in  northeast 
Colorado  that  was  well  observed  and  measured  during  the  1976  field  pro¬ 
gram  of  the  National  Hail  Research  Experiment  (NHRE) .  The  case  was 
chosen  primarily  because  of  the  availability  of  detailed  hydrometeor 
habit  and  size  data.  The  analysis  to  be  described  required  a  data  set 


ghted 


of  simultaneous  water  contents  and  reflectivities  along  a  flight  track. 


Ideally,  surface-based  radar  reflectivities  should  have  been  used  to 
match  with  simultaneous  aircraft  measurements.  However,  the  available 
radar  data  for  the  NHRE  case  were  insufficient  in  spatial  and  temporal 
resolution  to  allow  meaningful  comparisons.  The  reflectivities  were, 
therefore,  derived  from  the  same  airborne  hydrometeor  size  and  concen¬ 
tration  measurements  used  to  derive  the  hydrometeor  masses.  Thus,  the 
identical  resolution  between  the  M  and  Z  values  was  achieved. 

The  selected  hailstorm  occurred  on  25  July  1976  end  is  described  by 
Sanborn  (1979).  The  aircraft  data  were  collected  by  the  South  Dakota 
School  of  Nines  and  Technology  armored  T-28  aircraft  during  three  pene¬ 
trations  of  the  storm  complex.  Figure  3.7a-c,  reproduced  from  an 
analysis  of  the  storm  by  Hunter  (1980),  illustrate  the  flight  tracks 
through  the  storm  complex.  The  tracks  are  superimposed  on  representa¬ 
tive  radar  reflectivity  PPI  displays.  The  flight  level  averaged  about 
6.5  km  MSL,  with  the  temperature  between  -10  and  -15  C.  The  three  pene¬ 
trations  sampled  a  wide  variety  of  conditions. 

Hydrometeor  data  were  collected  with  two  particle  measuring  sensors 
on  board  the  T-28  aircraft:  a  Particle  Measuring  Systems  (PMS)  2-D 
cloud  probe  and  a  hail  spectrometer.  These  instruments  and  the  methods 
employed  to  reduce  the  data  at  NCAR  are  described  by  Heymsfield  and 
Parrish  (1979).  Table  3.2  gives  the  size  ranges  of  the  spectral  chan¬ 
nels  into  which  the  measured  particles  were  grouped.  Each  particle’s 
habit  was  classified  into  one  of  the  six  categories  listed  in  Table  3.3. 
For  particle  sizes  between  400  and  10,000  pm,  the  habit  was  determined 
objectively  from  the  2-D  images  and  microphysical  measurements,  as  out¬ 
lined  by  Heymsfield  and  Parrish  (1979).  Smaller  particles  were  assumed 
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Figure  3.7a:  Radar  reflectivity  at  1800:18  MDT  and  3.3°  elevation  angle 
showing  corrected  T-28  track  for  penetration  1  from  1754 
to  1805  MDT.  The  track  is  cross-hatched  where  >  5  m  s-* 
updrafts  were  encountered,  and  positions  are  narked  every 
minute  and  labelled  every  five  minutes.  The  altitude 
range  of  the  penetration  is  given  at  the  left  side  of  the 
figure.  A  'H'  denotes  a  location  where  hail  was  collected 
at  the  surface.  Reflectivity  contours  are  given  by  the 
legend  at  the  upper  right  and  arcs  of  constant  elevation 
are  labelled  in  kilometers  (from  Hunter,  1980). 


25  JULY  1976  T-28  PENETRATION  2 


Figure  3.7b:  Same  as  Fig.  3.7a  except  PPI  is  at  1816:07  MDT  and  track 
“j”  Penetration  2  from  1810  to  1820  MDT  (from  Hunter, 


Spectral  Size  Intervals  and  Assumed  Particle  Dimensions  from 
the  Concatenation  of  PMS  2-D  Cloud  Probe  and  Hail  Spectrometer  Data 
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to  be  plates  and  larger  particles  were  assumed  to  be  graupel/hai 1 .  Note 
in  Table  3.3  that  there  is  no  category  for  liquid  particles.  Recent 
research  with  PMS  2-D  data  by  Heymsfield  (1982,  personal  communication) 
on  NHRE  and  Oklahoma  cases,  and  our  own  experience  with  a  cumulonimbus 
case  in  Montana,  have  indicated  that  even  at  levels  of  near  -10°C, 
supercooled  drops  appear  in  the  2-D  images  and  should  be  classified  as  a 
separate  particle  category.  However,  their  sizes  and  concentrations  are 
sufficiently  small  that  their  contribution  to  total  water  mass  and 
reflectivity  is  negligible  in  contrast  to  the 

TABLE  3.3 

Habits  Used  for  Classification  of  Hydrometeors 

ID  Code  Habit  Symbol 

1  Graupel/Hail  GR 

2  Rimed  Dendrite  RD 

3  Rimed  Aggregate  of  Dendrites  RA 

4  Plates  PL 

5  Umrimed  Dendrite  UD 

6  Unrimed  Aggregates  of  Dendrites  UA 

7  Shedded  Particle  (Rejected) 

graupel  contribution.  Consequently,  no  category  for  supercooled  rain 
appears  in  Table  3.3. 

The  processed  T-28  data  from  NCAR  consisted  of  particle  concentra¬ 
tions  as  a  function  of  habit  and  spectral  channel  for  10-s  samples  along 
the  flight  tracks  in  Fig.  3.7,  Thus,  for  the  flight  speed  of  near  100 
m/s,  each  sample  was  collected  over  a  distance  of  about  1  km.  Particle 
mass  concentrations  (M)  were  calculated  using  the  following  relation- 


ship , 
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M  =  I  I  N 

i  j 


.  .m,  . 


(3.13) 


where  i  represents  all  channel  sizes  (Table  3.2),  j  represents  all 
habits  (Table  3.3),  N.,  is  the  particle  concentration,  and  ra..  is  the 
mass  of  an  individual  particle.  As  described  by  Heymsfield  and  Parrish 
(1979),  particle  mass  is  related  to  its  dimension  L  through  power  func¬ 
tions  of  the  form  m_  =  aL^,  with  each  habit  having  a  unique  set  of 
empirical  constants  a  and  b.  The  radar  reflectivity  factors  (Z)  were 
calculated  using  the  relationship  from  Heymsfield  and  Parrish 

Z  =  r  Z  N.  .  m*  x  3.648  x  109  x  (0.19/0.93)  (3.14) 

6  ij  1J  lj 

-1  -3 

where  N..  is  in  units  of  1  and  m..  is  in  units  of  g  m  .  The  factor 
ij  ij 

9 

3.648  x  10  is  a  unit-conversion  factor  relating  the  square  of  the 
particle's  mass  to  the  6th  power  of  its  melted  diameter.  The  ratio 
(0.19/0.93)  converts  the  reflectivity  factor  for  the  ice  particles  (Z) 
to  an  'effective  reflectivity  factor'  (Z^)  that  assumes  a  dielectric 
constant  for  liquid  water.  (Hereafter,  'reflectivity'  will  be  used  to 
mean  ’effective  reflectivity  factor.’) 

The  mass  concentrations  calculated  from  the  airborne  measurements 
using  Cq.  (3.13)  (M^)  are  plotted  as  a  function  of  time  in  Fig.  3.8a  and 
8b.  The  total  mass  concentrations  are  plotted,  along  with  the  mass  con¬ 
centrations  for  graupel,  plates,  and  rimed  dendrites.  The  mass  concen¬ 
trations  and  sample-to-sample  variability  of  rimed  aggregates,  unrimed 
dendrites  and  unrimed  aggregates  were  similar  to  the  values  for  the 
rimed  dendrite  plot;  these  additional  habits  are  omitted  from  the  fig¬ 
ures  for  clarity.  It  can  be  seen  in  Figures  3.8a  and  3.8b,  that  the 
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Fignre  3.8a:  Maaa  concentration*  of  particles  from  airborne 

measurements  on  25  July  1976,  Penetration  1.  Included  are 
plots  for  total  mass  concentration  and  contributions  due 
to  graupel,  plates,  and  rimed  dendrites. 
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variation  of  along  the  flight  track  is  large,  with  successive  10-s 
samples  frequently  differing  by  a  factor  of  5,  and  that  graupel  parti¬ 
cles  contribute  most  to  the  total 

The  reflectivities  calculated  from  the  particle  measurements  (Z^) 
using  (3.14)  are  plotted  as  dfiZ  versus  time  in  Pig.  3.9a  and  3.9b.  Con¬ 
siderable  variability  is  evident  in  the  reflectivity  values  in  the  fig¬ 
ures.  Hail  occurrence  near  1800  and  1815  MDT  caused  the  highest  reflec¬ 
tivities.  The  contribution  of  graupel  to  total  reflectivities  is  so 
large  that  the  graupel  reflectivity  values  are  not  plotted,  since  they 
nearly  coincide  with  the  total  reflectivity  curve.  The  contribution  due 
to  the  other  particle  habits  are  very  small,  typified  by  the  curve 
included  for  rimed  dendrites.  The  contribution  of  plates  to  the  total 
is  comparable  to  all  habits  except  the  dominant  graupel  category 
(Fig.  3.8).  However,  the  effect  of  plates  on  reflectivities  is  negligi¬ 
ble,  implying  large  concentrations  of  small  particles.  The  average  par¬ 
ticle  spectra  (from  penetrations  1,  2  and  3)  for  each  habit,  gr  phed  in 
Fig.  3.10,  confirm  the  large  number  of  small  plates.  Fig.  3.10  also 
illustrates  the  dominance  of  graupel  at  larger  sizes. 

In  order  to  study  the  factors  which  cause  variability  in  M-Z  rela¬ 
tionships,  a  least-square,  best-fit  power  function  of  the  form 

M.  =  aZ*  (3.15) 

A  Ae 

6  3 

was  derived,  where  Z^  is  the  effective-reflectivity  (mm  /m  )  deter¬ 
mined  from  the  airborne  particle  measurements  using  "4).  Ms  is  the 
mass  concentration  (g/m  )  derived  from  (3.13)  and  a  and  b  are  constants 
of  the  functon.  Included  in  the  sample  were  all  pairs  of  and  Z^e  for 

which  10  dBZ^  <  (lOlogZ^)  <  65  dBZ^,  or  160  samples.  These  pairs  are 
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Figure  3.9a:  Radar  effective-reflectivity  factors  calculated  from  the 
airborne  particle  measurements  from  25  July  1976, 
Penetration  1.  Included  are  plots  for  total  reflectivity 
and  contributions  due  to  rimed  dendrites.  The  plot  for 
graupel  is  omitted,  since  it  nearly  coincides  with  the 
total  plot. 
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Figure  3.10:  Average  particle  spectra  for  the  six  identified  habits 

from  Penetrations  1.  2,  and  3  on  25  July  1976.  Included 
in  the  spectra  are  the  160  samples  for  which  10  <  101ogZe 
i  ^  dBZe .  The  contribution  due  to  particles  grouped  into 
each  spectral  channel  is  plotted  at  the  mid-point  of  that 
channels’  spectral  range,  noted  at  the  top  (see  also  Table 
3.2) . 


plotted  in  Fig.  3.11  as  a  scatter-gran  of  logM  vs.  logZ^.  Solving  for 
the  best-fit  line  to  these  data,  a  and  b  in  (3.15)  become  0.065  and 
0.196,  respectively,  yielding  a  correlation  coefficent  of  0.50  and  a 
standard  error  factor  of  0.349.  This  best-fit  line  is  plotted  in  Fig. 

3.11. 

Now,  (3.15)  becomes 

Mn  =  0.065  Z°.‘196  (3.16) 

D  Ae 

where  is  a  derived  mass  concentration  that  lies  along  the  solid  line 
in  Fig.  3.11.  Values  of  Mp  were  calculated  using  this  expression  for 
each  value  in  Fig.  3.9  and  plotted,  along  with  the  corresponding 
value,  versus  time  in  Figs.  3.12a  and  3.12b. 

Coherent  periods  of  generally  under-estimated  M„  values  (M„  <  M  ) 

U  U  A 

occur  in  Figures  3.12a  and  3.12b  between  1755  and  1800  MDT  and  between 
1817  and  1823  MDT.  Slightly  less  coherent  periods  of  over-estimated 
values  (Mp  >  M^)  occur  between  1812  and  1817  MDT  and  after  1823  MDT. 

The  coherent  patterns  of  over-estimated  and  under-estimated  M  values 
suggest  that  the  variations  are  real  and  not  instrument  noise.  The  160 
(ma.  Z^  )  pairs  were  then  divided  into  two  groups,  Mp  >  M^  (below  the 
best-fit  Mp  line  in  Fig.  3.11)  and  Mp  <  (above  the  line  in  Fig. 

3.11),  in  order  to  determine  the  cause(s)  for  the  two  groups. 

T>  average  M^  and  Z ^  values  for  each  group  were  computed  for  each 
particle  habit  and  for  all  habits  combined.  The  results  are  graphed  in 
Fig.  3.13.  Values  of  M^  and  Z ^  for  all  habits  combined  are  plotted  in 
the  center  column  of  the  left  and  right  panels  in  Figure  3.12.  The 
values  of  M^  and  Z ^  for  the  over-estimated  (Mp  >  M^)  and  under¬ 
estimated  (Mp  <  M4)  groups  are  plotted  to  the  left  and  right, 
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lues  of  mass  concentrations  calculated  from  particle 
ectra  (*A) .  and  calculated  from  MD  vs.  ZZe  relationship 
D; *  Plotted  versus  time  for  Penetration  1.  Darkly 
aded  regions  between  the  plots  show  portions  of  the 
ight  track  where  M„  <  MA,  and  lightly  shaded  regions 
e  where  Mq  >  M». 
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Figure  3.12b:  Sane  as  Fig 


.  3.12a,  except  for  Penetrations  2  and  3 


REFLECTIVITY 


CONCENTRATION 


13:  Values  of  MA  (1  eft  panel)  and  ZAe  (right  panel)  for  the 
total  set  of  160  samples  (center  columns  of  each  panel), 
for  the  76  samples  where  Mp  >  MA  (left  columns),  and  for 
the  84  samples  where  Mp  <  MA  (right  columns).  Plotted  a 
values  for  each  habit  specified  in  Fig.  3.10  and  for  all 
habits  combined  ('T'). 
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respectively,  of  the  center  columns  in  Figure  3.13.  It  can  be  seen  in 
Figure  3.13,  amplifying  the  results  indicated  in  Figs.  3.8  and  3.9,  that 
graupel  (GR)  contributes  most  to  total  M^  and  Z ^  values  (T) .  The  con¬ 
tributions  by  the  other  five  habits  are  much  less,  1  to  2  orders  less 
for  and  4  to  6  orders  less  for  Z^ .  Plates  (PL),  due  to  their  large 
concentrations  of  small  sizes,  contribute  the  most  of  these  five  habits 
to  M^,  but  the  least  to  Z^.  The  dominance  of  graupel  to  and  2^ 
values  for  both  groups,  along  with  no  significant  variation  in  the  con¬ 
tribution  from  the  other  habits,  indicates  that  variations  in  the  grau¬ 
pel  size  distribution  must  account  for  the  systematic  differences 
between  the  M  >  M  and  M  <  M  groups. 

In  Fig.  3.14,  the  average  graupel  spectra  are  plotted  for  both  the 
Mp  >  and  Mp  <  groups  (the  average  graupel  spectrum  for  the  entire 
sample  in  Fig.  3.10  lies  between  the  spectra  here).  Also  plotted  in 
Fig.  3.14  are  the  spectral  M^  #nd  Z^e  values  for  the  groups. 

It  can  be  seen  in  Figure  3.14,  that  there  are  more  large  graupel 
particles  (d  >_  ~  6000pm)  in  the  Mp  >  group  than  in  the  Mp  <  group. 
Conversely,  there  are  fewer  small  graupel  particles  (d  <  6000pm)  in  the 
Mp  >  group  than  in  the  Mp  <  group.  These  particle  concentration 
variations  cause  an  increase  in  the  mass  distribution  for  particles  with 

d  >  ~  6000pm  and  a  decrease  for  particles  with  d  <  6000pm  in  the  M  >  M. 

D  A 

group,  relative  to  the  Mp  <  M^  group.  The  increase  in  mass  is  less  than 
the  decrease,  causing  a  net  decrease  in  the  total  mass  concentration 
<Ma); 

M  <  M 

<W  W 


as  seen  in  Figure  3.13. 


1! 
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Further,  because  of  the  variation  in  particle  concentration,  there 


is  a  corresponding  variation  in  reflectivity  values  between  the 

and  Mp  >  groups  as  shown  in  Figure  3.13:  the  group  shows  an 

increase  for  particles  with  d  >  6000  pm  and  a  decrease  for  smaller  par¬ 


ticles,  relative  to  the  group.  The  increase  in  reflectivity 

values  for  particles  with  d  >  6000pm  is  orders  of  magnitude  greater  than 
the  decrease  for  smaller  particles,  due  to  Z  being  proportional  to  d^ . 


Thus,  a  net  increase  in  Z  results  as  shown  in  Figure  3.13; 


Z.  >  Z 

<Vma>  (Vma> 

The  larger  Z  value  for  the  group  predicts  a  larger 

value  using  (3.16)  than  for  the  Mp  <  group.  However,  the  measured 
mass  concentration  was  greater  for  the  Mp  <  group  (Fig.  3.13).  The 
reason  for  this  discrepancy  can  be  seen  from  Fig.  3.14.  Given  the 
approximate  linear  nature  of  the  log-log  particle  concentration  distri¬ 
bution,  maximum  sensitivity  of  to  particle  concentration  variation 
occurs  at  sizes  below  6000pm,  peaking  near  3000pm.  On  the  other  hand, 

maximum  sensitivity  of  Z  occurs  at  the  largest  sizes  (>  6000pm). 

Ae 

Thus,  samples  with  relatively  large  numbers  of  particles  >  6000pm  and/or 
small  numbers  of  particles  near  3000pm  (represented  by  the  Mp  >  dis¬ 
tribution  in  Fig.  3.14)  have  larger  Z Ae  values  and/or  smaller  values 
relative  to  the  best  fit  Z-M  relationship,  giving  the  cluster  of  points 
below  the  line  in  Fig.  3.11.  Samples  with  relatively  large  numbers  of 
particles  near  3000pm  and/or  small  numbers  of  particles  >  6000pm 


(represented  by  the  Mp  <  M^  distribution  in  Fig.  3.14)  have  larger  M^ 
values  and  smaller  Z ^  values  relative  to  the  best-fit  relationship, 
giving  the  cluster  of  points  above  the  line  in  Fig.  3.11. 


/,q 

It  also  can  be  implied  from  an  examination  of  Figs.  3.11  and  3.14 

5  6  3 

that  the  samples  with  the  highest  reflectivities  (>  10  mm  /m  )  are  most 
likely  due  to  the  presence  of  hailstones  _>  1  cm.  Consequently,  if  the 
ten  or  so  largest  Z^  samples  are  eliminated  from  the  scattergram  in 
Fig.  3.11,  the  best-fit  line  to  the  remaining  data  wonld  lave  a  steeper 
slope,  0.256  instead  of  0.196. 

This  analysis  has  shown  that  the  graupel  size  distribution  was  most 
critical  to  the  variability  of  the  Z-M  relationship  in  a  northeast 
Colorado  hailstorm.  This  variability  is  due  to  different  spectral  sizes 
at  which  M  and  Z  variations  are  most  sensitive  to  change  in  number  con¬ 
centrations.  This  conclusion  supports  an  earlier  result  from  our  meas¬ 
urement  and  modeling  studies  of  Florida  cumulus.  Consequently,  accurate 
representation  of  graupel  formation  and  the  distribution  of  graupel  par¬ 
ticles  in  the  3-D  numerical  simulation  model  is  critical  to  proper 
reflectivity  predictions. 

The  total  mass  concentration  values  in  Figures  3.8a  and  3.8b  were 
utilized  along  with  the  corresponding  total  reflectivity  values  in  Fig¬ 
ures  3.9a  and  3.9b  to  compute  values  of  lOlogZ'/mass  concentration  for 
comparison  with  model-predicted  values  in  Figure  3.6.  The  computed 
ratios  are  illustrated  in  Figure  3.15  Two  results  are  apparent  upon 
comparing  the  ratios  in  Figures  3,6  and  3.15.  First,  the  variability  in 
the  model-predicted  and  measured  ratios  are  low.  Second  the  regions 
with  the  largest  ratios  are  at  cloud  edge.  These  results,  although 
based  on  a  limited  data  set,  are  in  qualitative  agreement.  This  agree¬ 
ment  indicates  that  the  predicted  water  content  and  reflectivity  values 
are  reasonable.  This  conclusion  is  consistent  with  those  made  from  our 
earlier  measurement  and  modeling  studies  of  Florida  cumulus. 


In  summary,  detailed  airborne  hydrometeor  measurements  were 
analyzed  for  a  northeast  Colorado  hailstorm.  The  analyses  focnsed  on 
determining  the  variations  in  the  relationship  between  water  contents 
and  radar  reflectivities  and  the  canse(s)  for  the  variations.  It  was 
found  that  systematic  variations  existed  in  the  hailstorm;  regions  where 
water-contents  calculated  from  radar  reflectivity  values  were  greater 
than  the  measured  contents  and  vice  versa.  The  variations  were  cansed 
by  relatively  greater  concentrations  of  the  largest  graupel  particles  (> 
6  ana  dia.)  combined  with  smaller  concentrations  of  particles  near  3  mm 
in  regions  with  over-estimated  water  contents,  and  lower  concentrations 
of  particles  larger  than  6  mm  combined  with  larger  concentrations  of 
particles  near  3  aim  in  regions  with  under-estimated  water  contents. 

The  hydrometeor  measurements  were  used  to  estimate  ratios  of  radar 
reflectivity  to  water  content  for  comparison  with  ratios  calculated  from 
3-D  model  simulations.  It  was  found  that  the  measured  and  calculated 
ratios  were  in  qualitative  agreement:  regions  of  high  ratios  occurred 
primarily  at  cloud  edges  where  a  few  large  graupel  particles  existed. 

The  results  of  this  study  demonstrate  the  importance  of  graupel 
particles  to  radar  reflectivity  measurement  and  calculations.  Conse¬ 
quently.  accurate  representation  of  graupel  formation  and  particle  spec¬ 
tra  in  the  3-D  numerical  simulation  model  is  critical  for  proper  reflec¬ 
tivity  predictions. 

3 .5  Explanations  of  Variations  in  M-Z  Relationships 

We  reported  in  section  3.4  that  systematic  variations  in  M-Z  rela¬ 
tionships  were  identified  that  were  caused  primarily  by  variations  in 
the  concentrations  of  graupel  particles  with  diameters  greater  than  6  mm 


and  with  diameters  near  3  mm.  We  found  that  variations  in  Z  were  highly 
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sensitive  to  variations  in  the  number-concentrations  of  the  graupel  par¬ 
ticles  2  6  sun,  while  variations  in  M  were  most  sensitive  to  variations 
in  the  number-concentrations  of  the  particles  near  3  mm.  We  reported  an 

_3 

average  relationship  (Eq.  3.16)  from  our  analysis  of  M  (g  m  )  =  0.065 

„0 . 1 96  .  6  -3, 

e 

As  a  result  of  these  findings.  Plank  (1982,  personal  communication) 
questioned  why  we  found  an  exponent  of  ~  0.2  instead  of  the  0.5  exponent 
used  in  our  3-D  model.  The  purpose  of  this  section  is  to  explore  that 
question  and  related  questions. 

It  can  be  seen  in  the  derivation  of  graupel  reflectivity  Z^  in 

Appendix  A  that  Z  a  M^,  where  the  graupel  mass  concentration  M  is 
8  Ntg  8  8 

a  product  of  air  density  and  the  graupel  mixing  ratio  r  .Values  of  r 

8  8 

are  the  basic  space-  and  time-dependent  variable  predicted  by  the  model 
and  drive  the  graupel  reflectivity  parameterization.  In  this  parameter¬ 
ization,  a  Marshall-Palmer  size  distribution  is  assumed  in  which  N  ,  a 

t8 

"total  graupel  number  concentration,"  must  be  specified.  In  the  South 

Florida  cumulus  simulations  in  section  3.3,  N  was  set  constant,  so 

tg 

that  the  Z-M  relation  indeed  had  an  exponent  of  0.5,  or  M  a  Z^*^. 

8  8 

Because  this  resulted  in  an  apparent  over-estimation  of  N  ,  especially 

tg 

at  early  stages  of  graupel  formation,  the  N  formulation  in  the  model 

*8 

was  revised  such  that  a  M^.  Thus,  the  model  M-Z  relation  became 
Mg  a  2*'°»  Jilh  >p  exponent  of  1.0  built  in  rather  than  0.5,  (Battan 
(1973)  cites  several  M-Z  relations  for  hail  in  which  the  exponent  is 


near  1.0.) 


In  the  following  analyses,  the  behavior  of  the  model  M-Z  relation¬ 
ship  (with  its  1.0  exponent)  is  compared  to  the  empirical  relationship 
derived  in  the  last  section  (with  its  ~  0.2  exponent).  The  M  and  Z 


values  from  the  sversse  measured  particle  concentration  spectrum  are 
compared  with  the  M  and  Z  values  from  a  model-parameterized  Marshall- 
Palmer  spectrum. 

We  consider  only  the  total-sample  average  spectra  in  the  analyses, 
since  the  differences  between  the  model  formulation  and  the  measured 
spectra  are  much  greater  than  those  between  the  different  measured  sam¬ 
ples.  While  not  explicitly  depicted  in  Fig.  3.14.  the  total-sample 
average  spectra  for  the  measured  particle  number  and  mass  concentrations 


and  reflectivities  all  lie  between  the  respective  and  Mp  <  M 

spectra. 


3.5.1 


The  graupel  particle  spectrum  in  the  3-D  model  is  after  the 


Marshall-Palmer  (M-P)  distribution: 


(3.17) 


-4 

where  N(D)  [cm  ]  is  the  particle  concentration  per  unit  size  interval 

_3 

at  particle  diameter  D  [cm],  N  [cm  ]  is  the  total  graupel  particle 

concentration,  and  D  [cm]  is  a  characteristic  diameter  of  the  distribu- 

o 

tion.  D  is  a  constant  of  the  distribution  which  must  be  specified, 
o 

For  the  Florida  reflectivity  studies,  the  constant  is  set  to  an 
equivalent  melted  diameter  of  0.054  cm  (related  to  a  characteristic 
dimension  of  a  M-P  raindrop  distribution  which  tends  to  stabilize  due  to 


drop  breakup) .  N  is  calculated  as  a  function  of  Dq  and  the  total 

3 

graupel’ mass  concentration,  M  [g/cm  ]: 
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(3.18) 


where  N  is  a  product  of  air  density  p  (g/cm  )  and  granpel  nixing  ratio 

g  a 

3 

r  (g/g),  p,  is  liquid  water  density  (1.0  g/cm  ).  and  0  is  the  charac- 
g  l  o 

teristic  equivalent  melted  particle  diameter  (cm).  Thus  N  is  a 

derived  total  particle  concentration,  such  that  if  all  N  particles  had 

a  melted  diameter  D  ,  the  total  mass  concentration  would  be  M  . 

o  g 

Because  radar  reflectivity  is  proportional  to  the  6th  power  of  an 
ice  particle's  equivalent  melted  diameter  (i.e.,  the  square  of  its 
mass),  the  formulation  for  reflectivity  in  the  model,  based  on  Eq. 
(3.17),  would  be  more  accurate  if  D  is  given  in  terms  of  melted  diameter 
(which  is  currently  not  done).  Thus,  modelled  spectral  values  of  grau- 
pel  reflectivity,  Z^,  as  well  as  spectral  mass  concentrations,  IT, 
depend  on  the  assumed  values  of  ice  density,  p^,  at  the  graupel  size  . 
However,  as  will  be  shown,  Zj  and  Mj  are  relatively  insensitive  to 
differences  in  liquid  vs.  ice  water  density,  compared  to  more  dominant 
assumptions  concerning  Ntg  and  D0. 

To  examine  how  the  model  would  simulate  the  25  July  1976  storm 
analyzed  in  section  3.4.2,  we  solve  for  N  in  Eq.  (3.18)  using  the 

3 

measured  average  total  mass  concentration  H  =  0.25  g/m  and  D  =  0.027 

g  o 

-2  -3 

cm  to  get  N  =  2.43  x  10  cm  Using  these  N  and  D  values  and  the 
tg  *  tg  o 

mid-channel  D  values  defined  by  the  various  channels  of  the  2-D  cloud 
probe  and  hail  spectrometer  reported  in  Table  3.2,  N(D)  values  were  cal¬ 
culated  using  Eq.  (3.17).  The  results  are  depicted  in  Fig.  3.16.  The 
result  using  a  frozen  diameter  directly  in  Eq.  (3.17)  is  seen  in  Fig. 
3.16  as  a  straight  line  on  the  log-linear  plot,  a  characteristic  of  a 
M-P  distribution.  The  other  two  lower  curves  of  the  three  lower  curves 
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Figure  3.16:  Plot  of  various  Marshall-Palmer  formulations  of  graupel 
size  distributions.  Results  of  calculations  using  Eq. 
(3.17)  as  in  the  current  3-D  model  simulations,  and  an 
alternate  form  of  Eq  (3.17). 
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in  Fig.  3.16  are  derived  from  Eq.  (3.17)  using  equivalent  Belted  diame- 
ters  corresponding  to  the  abscissa  values  of  frozen  spectral  diameter  D. 

3 

The  plot  for  melted  diameters  using  p^  =  0.6  g/cm  ,  as  assumed  by  the 

model,  is  a  straight  line,  since  the  constant  density  affects  D  by  a 

constant  factor.  The  plot  for  melted  diameters  using  the  variable  p  is 

g 

based  on  the  frozen-to-melted  diameter  assumptions  described  by  Heyms- 

field  and  Parrish  (1979).  Their  iormulation  of  density  ranges  from  less 

3  3 

than  0.1  g/cm  for  the  smallest  graupel  particles  to  just  over  0.6  g/cm 

for  the  largest.  This  curve  would  also  be  a  straight  line  if  plotted 

against  the  equivalent  melted  diameters. 

The  number-concentration  of  particles  in  a  discrete  size  interval 

(or  sensor  spectral  channel),  is  given  by  the  integral  of  the  curve 

in  Fig.  3.16  over  that  size  interval.  Integrating  Eq.  (3.17), 
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(3.19) 


The  three  lower  curves  of  N(D)  in  Fig.  3.16,  which  were  based  on  the 

current  model  formulation  of  N  ,  were  integrated  over  the  spectral 

tg 

channel  limits  D^  and  D^  specified  in  Table  3.2,  using  the  actual  frozen 
diameter  or  the  appropriate  equivalent  melted  diameter  as  required  for 
each  curve.  The  resultant  'current  model  formulation'  particle  concen¬ 
tration  spectra  are  plotted  in  Fig.  3.17,  along  with  the  'aircraft- 
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Figure  3.17: 


Particle  Diameter  tyxm) 

Plot  of  particle  onuber  concentration*  integrated  over 
spectral  aise  channels.  Measured  hydrometeor  size  spectra 
■  .  ;  current  3-D  model  M-P  spectra  from  Eq .  (3.19)  using: 

frozen  diameters  ,*  ,  melted  diameter*  (variable  density), 
x  ,  melted  diameters  (pg  =  0.6  g  cm~3),  4  ;  alternate  M-P 
formulation  using:  melted  diameters  (variable  Pg).S.  and 
melted  diameters  (o.  “  0.6  «  cm~3) .  yf  . 
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measured'  ('real')  spectrum.  It  can  be  seen  in  Fig.  3.17  that  the 
differences  due  to  integrating  between  discrete  "frozen*'  size  inter¬ 
vals  vs.  the  two  melted  size  intervals  (nsing  two  different  ice  density 
assumptions)  are  much  less  than  the  difference  between  any  of  the  three 
curves  and  the  measured  spectrum.  Further,  smaller-sized  particles  are 
over-estimated  in  number  by  about  an  order  of  magnitude  in  comparison 
with  the  measured  spectrum  and  the  mid  thru  largest  sizes  are  very  seri¬ 
ously  under-estimated.  Due  to  the  large  magnitudes  of  the  over¬ 
estimated  concentrations  of  small  particles,  the  modelled  total  particle 
concentrations  are  several  times  larger  than  the  measured  total  concen¬ 
tration  (see  total  concentrations  summarized  at  the  bottom  of  Fig. 

3.17) . 

3.5.2  Model  Formulation  of  Graunel  Mass  Concentration  and 
Reflectivity  Spectra 

The  modelled  mass  concentration  spectra,  corresponding  to  the 
'current  model  formulation'  particle  concentration  spectra  in  Fig.  3.17 
are  depicted  in  Fig.  3.18,  along  with  the  aircraft-measured  mass  concen¬ 
tration  spectrum.  The  appropriate  ice  density  assumptions  as  assumed  by 

Heymsfield  and  Parrish  (1979),  variable  p  )  and  as  used  in  the  model  (p 

8  8 
3 

=  0.6  g/cm  ),  respectively,  were  applied  to  produce  the  two  ’current 
model  formulation*  curves  (the  'frozen  diameter'  curve  corresponding  to 
the  third  'current  model  formulation’  curve  in  Fig.  3.17  has  been  omit¬ 
ted,  since  using  either  p  assumption  at  this  point  produces  no  signifi- 

8 

cant  differences  from  the  other  curves).  As  can  be  seen  in  Fig.  3.18, 
differences  between  the  modelled  and  measured  mass  concentration  spectra 
are  similar  to  the  differences  in  the  particle  concentration  spectra  in 
Fig.  3.17:  mass  in  the  smaller  size  channels  is  over-estimated  with 


respect  to  the  aircraft  measurements,  while  at  the  large  sizes,  mass  is 
very  much  under-estimated.  Due  to  the  large  magnitudes  of  the  over¬ 
estimates  at  the  smaller  sizes,  the  total  mass  concentrations  are  over¬ 


estimated  by  several  times  (see  total  mass  concentration  summaries  at 
the  bottom  of  Fig.  3.18). 

The  effect  of  the  differences  between  modelled  vs.  measured  parti¬ 
cle  concentration  spectra  shown  in  Fig.  3.17  on  the  'cnrrent  model  for¬ 
mulation'  vs.  measured  radar  reflectivity  spectra  is  shown  in  Fig.  3.19. 
It  can  be  seen,  that  the  reflectivities  calculated  from  the  'current 
model  formulation’  particle  concentration  spectra,  using  the  two 
assumptions,  vary  insignificantly  from  each  other,  compared  to  the  vari¬ 
ation  between  either  modelled  spectrum  and  the  aircraft-measured  reflec¬ 
tivity  spectrum.  The  modelled  reflectivity  spectra  have  similar  trends, 
relative  to  the  measured  spectrum,  as  are  exhibited  by  the  modelled  par¬ 
ticle  and  mass  concentration  spectra:  the  spectral  Z^  values  are  over¬ 
estimated  at  the  small  particle  sizes  because  of  the  over-estimated 
number  of  small  particles,  while  the  reverse  is  true  at  the  large  parti¬ 
cle  sizes.  However,  the  net  effect  on  the  modelled  total  reflectivity 
is  opposite  from  the  net  effect  on  the  modelled  total  mass  concentra¬ 
tion.  Since  the  dependency  of  Z(D)  causes  Z  to  be  highly  sensitive 
to  large  particles,  the  over-estimate  of  Z  with  respect  to  the  measured 
values  at  small  sizes  (due  to  the  over-estimate  of  particle  concentra¬ 
tions  at  small  sizes)  is  more  than  offset  by  the  under-estimate  of  Z  at 
large  sizes  (due  to  the  many-orders  under-estimate  of  large-particle 
concentrations) .  Indeed,  the  aircraft  Z,  values  peak  at  the  large 
sizes,  where  the  corresponding  modelled  Z^  values  are  negligible.  The 
resultant  total  modelled  reflectivity  is  an  order  of  magnitude  less  than 


Reflectivity  t  d  BZe ) 


Figure  3.19:  Measured  and  modelled  spectral  radar  reflectivities 
Symbols  are  the  same  as  in  Fig.  3.17. 
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that  measured  (see  total  reflectivity  summary  at  the  bottom  of  Fig. 
3.19)  . 


Thus,  the  model  formulation  of  a  M-P  distribution  using  the  con¬ 
stants  as  specified  in  Eq.  (3.18)  result,  for  this  comparison  with  meas¬ 
ured  hydrometeor  spectra,  in  an  over-estimate  of  total  graupel  mass  con¬ 
centration  and  an  under-estimate  of  total  graupel  reflectivity,  relative 
to  the  measured.  Hence,  either  of  the  modelled  parameters  a  and/or  b 

(where  b  =  1.0,  since  N  a  II  )  in  the  modelled  M  =  aZ*5  relation  would 

tg  g 

appear  to  be  too  large. 

3.5.3  Alternate  Model  Formulation  of  Graupel  Spectra 
What  is  clearly  needed,  as  seen  in  Fig.  3.17,  is  a  modelled  distri¬ 
bution  with  fewer  small  particles  and  more  large  ones.  The  question  is, 
can  a  M-P  distribution  do  this?  A  more  representative  distribution  can 
be  accomplished,  to  some  extent,  by  modelling  a  distribution  with  a 
flatter  slope  and  a  smaller  intercept  in  Fig.  3.16;  i.e.,  from  Eq. 

(3.17),  a  larger  value  of  D  and/or  a  smaller  value  of  N  .  It  should 

o  tg 

be  noted  that  the  previous  value  of  N  was  based  on  studies  of  raindrop 

tg 

spectra,  and  that  this  variable  has  yet  to  be  ’calibrated'  for  graupel 

spectra.  Thus,  to  accomplish  the  desired  adjustments,  we  based  a  simple 

alternate  formulation  of  N  on  the  measured  total  particle  concentra- 

tg 


tion  (over  all  spectral  channels)  from  Fig.  3.17.  The  new  value  of  N 


tg 


-3  -3 

is  2.05  x  10  cm  The  value  is  an  order  of  magnitude  less  than  the 

value  from  Eq.  (3.18)  which  was  based  on  the  measured  mass  concentra¬ 
tion.  Instead  of  assuming  a  characteristic  particle  size  D  of  0.027 

o 

cm,  as  before,  we  now  compute  from  the  measured  total  particle  and 


mass  concentrations: 


v»;.^  .  1J22 , 

1  ON  -  A~3  -3 

tg  2.05  x  10  cd 
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D  =  0.0616cm  . 

o 


(3.20) 


Note  that  we  have  reformulated  N  and  D  in  terms  of  an  average 

tg  o 

measured  particle  concentration  spectrum.  Since  the  model  predicts  only 

r  as  a  basic  variable  and  not  the  particle  concentration,  this  reformu- 
8 

lation  would  not  be  able  to  be  implemented  into  the  model.  Furthermore, 

since  no  explicit  formulation  of  N  in  terms  of  M  is  provided  by  this 

tg  8 

1  2 

illustrative  reformulation,  its  effect  on  the  modelled  Z  a  ~ —  M  rela- 

g  N  g 

8  tg  b 

tion  (and  thus  on  the  exponent  in  the  M  =  aZ*5  relation)  cannot  be  exam¬ 
ined.  It  serves  here  only  to  modify  the  constants  N  and  D  (and  thus 

tg  o 

the  slope  and  intercept)  in  the  modelled  Marshall-Palmer  distribution. 
The  modification  yields  reformulated  spectral  M  and  Z  values  that  can  be 
compared  to  the  first  formulation. 

The  resulting  alternate  M-P  distributions  are  shown  in  Fig.  3.16 
(using  the  identical  ice  density  assumptions  used  for  the  various  curves 
in  the  first  formulation).  As  seen  in  Fig.  3.17,  integration  of  the  new 
N(D)  curves  over  discrete  intervals  clearly  gives  reduced  concentrations 
in  the  smaller  channels  and  larger  concentrations  at  the  larger  channels 
than  with  the  first  formulation.  The  resultant  net  concentrations  are 
much  closer  to  the  measurements  (see  bottom  of  Fig.  3.17).  The  measured 
curve  is  better  reproduced,  with  the  errors  now  being  under-estimation 
at  the  smallest  sizes,  over-estimation  at  mid  sizes,  and  reduced  under¬ 
estimation  of  large  particle  concentrations. 

The  spectral  mass  concentrations  resulting  from  the  alternate  M-P 


spectra  are  shown  in  Fig.  3.18.  It  can  be  seen  in  Fig.  *.18,  thst  the 
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region  of  over-estimation  has  been  shifted  to  the  middle  part  of  the 
spectrum  (still  an  order  of  magnitude  too  high),  and  the  M^  values  at 
the  large-particle  sizes  remain  nnder-estimated.  The  total  mass  concen¬ 
trations  are  corrected  slightly  in  the  right  direction,  but  the  values 
remain  several  times  larger  than  the  measured  values. 

The  spectral  reflectivities  (Fig.  3.18)  calculated  from  the  alter¬ 
nate  M-P  spectra  show  improvement.  The  total  Z  values  are  now  within  a 
factor  of  2  of  the  values  calculated  from  the  airborne  hydrometeor  spec¬ 
tra.  Note  that  the  increase  in  total  Z  has  been  achieved  by  over¬ 
estimating  Z^  values  resulting  from  mid-size  particles  between  0.2  to 
1.0  cm,  while  still  grossly  under-estimating  the  most  significant  meas¬ 
ured  Z^  values  resulting  from  particles  with  sizes  >  1  cm. 

The  alternate  M-P  formulation  results  in  increased  Z  and  decreased 
M  (both  towards  the  measured  values),  suggesting  a  reformulated  Z-M 
relationship  that  is  closer  in  line  to  the  measured  0.2  exponent.  It 
should  again  be  stressed  that  the  model  does  not  have  a  0.5  exponent 
'built  in',  because  while  Z  is  proportional  to  the  mass  concentration 


squared  (inverse  of  0.5),  it  is  also  inversely  proportional  to  N  , 
which  is  also  related  to  mass  concentration  either  directly  [as  in 


(3.18)]  or  indirectly  through  some  sort  of  D  parameterization. 

o 

For  this  particular  case,  it  seems  unlikely  that  the  measured  pai 


tide  number-concentration  'tail*  at  the  large  end  of  the  spectrum  in 


Fig.  3.17  can  be  modelled  with  a  M-P  distribution,  regardless  of  the 
adustmenta  to  Ntg  and  D0.  The  'tail'  is  due  to  the  presence  of  large 
hail  in  this  particular  storm,  where  a  raindrop-type  breakup  mechanism 
would  not  provide  a  stable  anchor  to  the  large-particle  end  of  the  grau- 
pel  spectrum.  Hence  the  hydrometeor  spectra  may  defy  M-P  description. 
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3.5.4  Cone  Ins  ion* 

The  M-Z  relationship  derived  in  the  section  3.4,  M  =  0.065Z®'*^, 

e 

resulted  in  an  exponent  that  is  significantly  less  than  the  0.5  to  1.0 
value  often  cited  in  the  literature.  The  0.2  exponent  was  dne  to  the 
presence  of  large  hail  in  some  of  the  samples. 

The  measured  hydrometeor  size  distribution  was  compared  with  the 
Marshall-Palmer  (M-P)  approximation  currently  formulated,  based  on  con¬ 
stants  N  and  D  ,  in  the  3-D  model.  It  was  found  that  these  constants 
tg  o 

have  much  more  significant  impact  on  M  and  Z  than  do  the  assumptions 
concerning  ice  density.  With  the  current  M-P  formulation,  it  was  found 
that  particles  <  2.5  mm  diameter  are  over-estimated  in  number  by  about 
an  order  of  magnitude,  in  comparison  with  the  measured  size  distribu¬ 
tion.  and  the  particles  >  2.5  mm  are  seriously  under-estimated.  These 
differences  produce  corresponding  over-estimations  and  under-estimations 
at  the  small-particle  and  large-particle  ends,  respectively,  of  both  the 
modelled  mass  concentration  and  reflectivity  spectra,  relative  to  the 
measured  M  and  Z  spectra.  However,  total  M  is  more  sensitive  to  the 
over-estimations  of  the  small-particle  concentrations,  while  total  Z  is 
most  sensitive  to  the  nnder-estimations  of  large-particle  concentra¬ 
tions.  Thus,  the  M-P  distribution,  using  N  and  D  constants  as 

tg  o 

presently  specified  in  the  3-D  model,  results  in  an  over-estimate  of 
total  mass  concentration  and  an  under-estimate  of  total  reflectivity, 
relative  to  the  measured  total  M  and  Z  values.  Hence,  an  exponent  of 
0.5  to  1.0  appears  to  be  too  large  for  the  M  and  Z  values  resulting  from 
analysis  of  northeast  Colorado  hailstorm  data. 

Accordingly,  an  illustrative  alternate  M-P  distribution  was  derived 


nsing  measured  total  particle  concentrations  for 


(which  cannot  be 
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explicitly  predicted  by  the  model),  and  Dq  was  computed  using  N  and 
measured  total  mass  concentration.  The  alternate  N-P  formulation 
resulted  in  increased  Z  values  and  decreased  M  values,  suggesting  a 
modelled  M-Z  relationship  closer  to  agreement  with  the  0.2  exponent 
resulting  from  the  hailstorm  analysis. 

For  this  particular  case,  it  seems  unlikely  that  the  measured  par¬ 
ticle  number-concentration  'tail'  at  hail  sizes  >  1.4  cm  can  be  modelled 

with  a  M-P  distribution,  regardless  of  the  adjustments  to  N  and  D  . 

tg  o 

Consequently,  improved  formulations  of  graupel  spectra  are  needed. 

M  vs.  Z  relationships  were  ivestigated  in  three  climatic  regions 
during  this  study.  In  section  3.2,  reflectivities  were  predicted  for  a 
South  Park,  Colorado  thunderstorm  using  existing  M-Z  relationships  in 
the  3-D  numerical  model.  In  section  3.3,  reflectivities  and  particle 
mass  concentrations  were  measured  in  Florida  thunderstorms  but  too 
coarsely  to  develop  rigorous  M  vs.  Z  relationships.  In  section  3.4  and 
3.5,  M-Z  relationships  were  derived  from  airborne  particle  size- 
distribution  measurements  and  a  size  distribution  -  Z  relationship  from 
Heymsfield  and  Parrish  (1979).  It  can  be  seen  from  these  investiga¬ 
tions,  no  consistent  analysis  was  applied  to  the  three  climatological 
regions.  Consequently,  it  was  not  possible  to  draw  conclusions  as  to 
the  geographic  or  climatic  contributions  to  the  variability  of  M-Z  rela¬ 
tionships.  What  can  be  said  from  these  analyses,  through,  is  that 
regions  where  low  number  concentrations  of  large  hail  stones  are  pro¬ 
duced,  anomalous  M-Z  variations  will  occur. 
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3.6  Appendix  A 


Radar  Reflects 


by  William  R.  Cotton 


In  our  analysis  of  radar  reflectivity  and  its  variability  with  liquid 
water  content,  we  shall  apply  the  Rayleigh  approximation  to  the  radar 
equation.  Thus  the  power  returned  to  the  radar  from  an  ensemble  of  cloud 
precipitation  elements  is 


P~  =  — 
r  g 


‘ E.,' 


The  terms  in  [  ]  are  all  properties  of  the  radar  and  distance  r 

from  the  radar  and  are  defined  in  Battan  (1973).  For  our  purposes,  the 

remaining  terms  are  the  main  area  of  investigation.  The  term  |K.|  is  a 

function  of  the  complex  index  of  refraction  of  the  scattering  elements. 

2 

For  an  all  water  cloud,  | K J  is  approximately  0.93  and  for  ice  particles, 

2 

| K |  is  0.19.  In  our  analysis  of  model  results,  we  shall  assume  that 
ice  particles  in  the  melting  zone  are  water-coated  and  therefore  have  a 
|K|2  of  0.93. 

Normally,  meteorologists  employ  the  so-called  reflectivity  factor  Z 
to  designate  the  effects  of  the  particle-size  distribution  on  the  reflected 
power  returned  to  the  radar, 


z  *  E  ni  Di6  • 

vol 

In  our  analysis,  we  shall  define  a  complex  index  of  refraction-weighted 
reflectivity  Z  such  that 

Z  =  |K| 2  2  n±  b±6  . 

vol 


1 


r 
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Reflectivity  factor  Z  for  Rain,  rr 
It  is  assumed  that 


N  / 

♦  (r)  -  -jr  e  m 
m 
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Thus,  Z  increases  linearly  with  the  rainwater  content  r  .  It  is 

r 

also  interesting  to  note  that  Z  has  a  cubic  dependence  on  the  assumed 

characteristic  radius  r  of  the  raindrop  distribution. 

m 


Reflectivity  factor  for  Graupel  r 


g 


For  graupel,  Z  is 


00 

•ik/J 


D6N(D)e"XDdD 


since 


N(D)  =  N  \e~XT)  , 
tg 


then 


Z  =  K 


i'2  v/D 


6  -*D,n 

e  dD  , 


or  letting  x  =  AD  =>  dx  =  XdD  , 


CO 

lKl|2  J  x6e“Kd>= 


lKil2V 


r(7) 


But 


1/  9  9 

tip  N  \  '3  ,  p  r 

x  =1-^-)  ? 

Parg  I  X  it  P  N 

'  g  tg 


Thus 


r(7)P. 


z  = 


K.  r 


g  "  2  2  ‘‘'i1  g  ’ 

6  it  p  N 

g  tg 


where  N.  =  p  r  /ti6  p  D 

tg  a  g  l  mg 
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Thus,  for  graupel ,  Z  is  proportional  to  and  inversely  proportional 


to  N 


Reflectivity  factor  for  Ice  Crystals  r^ 
For  ice  crystals,  Z  is  given  by 


Z  =  |K. |2  N.d.6 
l  11 


where  is  the  crystal  concentration  defined  externally  in  many  cases 
by  the  Fletcher  formula,  and 

L  —7 

d^  =  ,  where  k^  =  0.515  mgm  for  <  1.7  x  10  gm,  and 

k^  =  0.192  mgm  ^  for  1.7  x  10  ^  gm  <  m ,  <  1 x  10  ^ 

and 


0.417 


and 


d.  =  k,  m. 
l  1  1 


Va 
"i  N. 

l 


for  m.  >  1  x  10  5  gm,  and 
k^  =  8.89  x  10  ^  mgm 


Since  is  nearly  always  proportional  to  nu  ,  this  shows  that  for 


ice  crystals 


12 

Z.  a  r.  /N. 

i  l  i 


Summary 

In  summary,  depending  upon  the  particular  form  of  water  substance  and 
the  assumed  distribution  functions  for  these  respective  components,  the 
relationship  between  reflectivity  factor  Z  and  water  content  varies  as  follows 
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For  rain  only 


Z  a  |K  | 2  r  3  r  , 

r  '  l 1  m  r  * 


For  graupel  only 


Z  a  K. 
g  1  i 


tg 


For  ice  crystals 


Zi  * 


2  r.' 

i 


In  the  supercooled  portions  of  a  simulated  cloud,  we  see  that  Z  varies 

substantially  with  total  precipitation  water  rp(rp  =  +  r  +  r.)  depending 

on  the  relative  contributions  of  r  ,  r  and  r..  We  should  also  note  that 

r  g  i 

2 

| K . |  will  be  assumed  equal  to  that  of  pure  ice,  except  below  the  melting 

2  2 

level  where  we  assume  that  Ik.I*"  =  Ik  ! 

1  l 1  1  C. 1 

To  analyze  the  variability  of  Z  with  r  ,  we  shall  define  Z  as  the 

P  r 

reflectivity  factor  computed  as  if  all  precipitation  were  in  the  form  of 
rain, 


Z  a  |K J2  r  3  r 
r  '  l '  ra  p 


Also,  we  define  Z  =  Z  +  Z  +  Z . , 
p  r  g  i’ 


and  then  plot  in  our  analyses. 


1 

—£■  along  with  plots  of 

r 


Z  and  r  , 
P  r’ 


i’ 


r 

P 
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4.0  Three-Dimensional  Cloud/Mesoscale  Model  Investigations  of  Potential 
Aircraft  Icing  Regions 
4.1  Background 

The  problem  of  predicting  aircraft  icing  conditions  is  primarily 
one  of  predicting  the  amount  of  supercooled  liquid  water  content  in  a 
cloud.  Simple  as  this  problem  may  seem,  it  actually  represents  a  very 
challenging  meteorological  problem.  This  is  because  the  amount  of 
supercooled  liquid  water  content  is  the  residual  of  the  differences 
between  large  production  and  removal  mechanisms.  Liquid  water  is 
primarily  produced  by  adiabatic  cooling  and  occasionally  by  radiative 
cooling.  However,  liquid  water  is  removed  from  the  cloud  by  entrainment 
processes  and  precipitation  processes.  In  a  supercooled  cloud  this 
means  removal  by  supercooled  rain,  vapor  deposition  growth  of  ice 
crystals  (i.e.,  via  the  Bergeron-Pinde isen  mechanism),  accretion  or 
riming  growth  of  ice  particles  including  pristine  ice  crystals,  graupel 
particles  and  aggregates  of  snow  crystals,  and  the  precipitation  of 
those  particles. 

In  an  earlier  simulation  of  orographic  clouds.  Cotton  et  al .  (1982) 
concluded  that  an  ice-crystal  aggregation  model  was  needed  to  properly 
predict  observed  amounts  of  precipitation  and  supercooled  liquid  water. 
Moreover,  the  ice-phase  competition  for  supercooled  water  is  greatly 
dependent  upon  the  concentration  of  ice  crystals.  Thus,  two  new 
developments  in  the  microphysical  model  are  needed  to  better  predict  the 
conditions  suitable  for  aircraft  icing  conditions.  These  are: 
i)  an  ice  crystal  aggregation  scheme 
ii)  prediction  of  the  concentration  of  ice-crystals  due  to  primary 
and  secondary  nucleation,  and  removal  by  collection  by  other 
precipitation  species  and  by  precipitation. 


f 

L 


The  rate  of  collection  among  a  homogeneous  population  of  ice 


crystals  is  given  by 


dN. 
_ 1 

dt 


CN 


=  -  K.  NT 

1  1 


(4.1) 


where  represents  the  concentration  of  "pristine" 
crystals  in  onr  model  and  is  the  collection  cross 
The  conversion  rate  of  ice  crystal  mixing  ratio 


(non-aggregated) 

section. 

r.  to  aggregates 


CN 


ia 


po  dt. 


=  +  K.  N.  r. 

ill 


(4.2) 


where  the  average  crystal  mass  m.  is  given  by 


m . 
i 


l 


o 


(4.3) 


and  is  the  air  density.  [Note  in  our  sign  convention  CN  indexed  ia 
(CNi#)  represents  a  positive  contribution  to  aggregates,  whereas  CNai 
represents  a  negative  contribution.]  A  pure  sedimentation  model  for 
aggregation  applied  to  our  homogeneous  population  of  pristine  crystals 
gives  =  0.  Therefore,  to  estimate  we  adapt  Passarelli  and 
Srivastava's  (1978)  stochastic  collection  kernel  model  which  estimates 
based  on  a  distribution  of  particle  densities  for  equal-sized 
crystals.  This  model  reduces  to 


J 
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K 


i 


hd: 


V.  E.X 
1  1 


(4.4) 


where  X  is  a  measure  of  the  variance  in  particle  fall  speed.  D^  and  V. 
are  the  average  crystal  diameter  and  fall  speed,  respectively. 
Passarelli  and  Srivastava's  best  estimate  of  X  is  0.25,  which  we  will 
use  as  a  first  guess  on  z.  We  plan  to  experimentally  adjust  z  on 
specific  observed  cases  to  calibrate  K.. 

4.2.2  The  Distribution  of  Aggregates 
Based  on  data  reported  by  Rogers’  (1974)  and  personal  discussions 
with  Passarelli  we  assume  aggregates  are  distributed  in  the  form 


N(D) 


N_  -D/D 
l  m 

7T*  e 


(4.5) 


where  N(D)  represents  the  spectral  density  of  aggregates  of  diameter  D, 

N  is  the  total  aggregate  concentration  and  D  is  a  "characteristic" 
l  m 

diameter  of  the  aggregate  population.  From  Roger's  (1974)  data,  we 
estimate 


D  =  0.33  cm 
m 

Employing  the  aggregate  density  relationship  given  by  Passarelli 
and  Srivastava  (1978): 

p,.  ■  c,  “'°'s  <«•<> 

-2.4 

where  {3^  =  0.015  g  m  and  integrating  over  the  entire  distribution 
(4.5)  ,  we  find 

N_  =  0.641  pi1  p  D~2'4  r 
T  lorn  a 


(4.7) 
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where  represents  the  mixing  ratio  of  aggregates. 
The  terminal  velocity  of  aggregates  is  given  by 


V 

a 


1/2 

ill  r-l/2 

iPoJ  ° 


1/2 

Pi. 


,1/2 


(4.8) 


where  g  is  the  acceleration  dne  to  gravity,  and  is  the  drag 
coefficient  here  taken  to  be  C  =  1.3.  Substitution  of  (4.6)  into  (4.8) 
gives  us: 


V 

a 


[8  P, 


LpoCD. 


1/2 


,0.2 


(4.9) 


Also  snbstitution  of  (4.6)  into  (4.5)  gives  us 

-l  i  a  'D/D 

N(D)  =  0.641  p/  p  r  D  3  e  m  (4.10) 

loam 


4.2.3  Aggregates  Collecting  Ice  Crystals 
Once  "pristine"  ice  crystals  have  been  converted  to  aggregates 
distributed  according  to  (4.10),  the  aggregates  can  collect  "pristine" 
ice  crystals.  The  change  in  concentration  of  crystals  by  aggregation  is 
given  by 


-  K(D  ,D  )  N  N(D  )  d  D 

a  i  i  a  a 


(4.11) 


For  large  aggregates  collecting  pristine  ice  crystals,  we  can 
employ  the  standard  gravitational  collection  kernel. 
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K  (D  ,D.) 

a  1 


S 

4 


(D 


Dir 


|V  -  V.|  E(a/i) 
a  l 


(4.12) 


where  we  approximate 


f V  <D  )  -  V. j  ~  |V  -  V. f 

1  a  a  i'  a  i 


The  average  terminal  velocity  of  ice  aggregates  is  given  by 


V 

a 


/  V  (D  )  m(D  )  N(D  )  dD 

a  a  a  a  a 

o _ 

r  P 
a  ro 


(4.13) 


which  after  substitution  of  (4.8),  (4.10),  and  (4.6)  becomes: 


V 

a 


=  2.49  0 


0.5 

1 


(4.14) 


We  can  now  estimate  the  change  in  aggregate  mixing  ratio  due  to 
collection  as 


CL. 

la 


m.  dN. 
i  _ l 

p  dt 
o 


=  ~  /  K(D  ,D . )  N.  N(D  )  d  D 
CL  Ni  o  8  1  1  a 


(4.15) 


Substitution  of  (4.12),  (4.14),  and  (4.10)  into  (4.15)  and  integration 
we  find 

CL.  =  0.503  0,  p  E(a/i)  |v~  -  \T|  D_2'4x  (4.16) 

la  lo  aim 

(2  D2  +  2  D.  D  +  D2)  r  r. 
m  i  m  x  a  l 

We  have  experimented  with  several  models  of  collection  efficiency 
E(a/i)  and  E^  in  (4.16)  and  (4.4)  based  on  data  reported  by  Passarelli 
and  Srivastava  (1978),  (1974),  Hallgren  and  Hosier  (1960)  and  Hosier  and 
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Hallgen  (1960).  The  temperature-influences  on  collection  efficiency  are 
evaluated  from  diagnosed  surface  temperatures  for  crystals  growing  by 
riming  and  vapor  deposition.  We  have  also  assumed  the  hydrodynamic 
collision  efficiency  is  1.4  as  evaluated  in  one  case  by  Passarelli  and 
Srivastava  (1978)  over  the  temperature  range  -12  to  -15°C.  Sensitivity 
experiments  with  these  models  demonstrated  that  the  Hosier  and  Hallgren 
model  results  in  an  underestimate  of  aggregation  formation  and  an 
incorrect  vertical  distribution  of  aggregates  where  E  ~  1.4  results  in  a 
realistic  vertical  distribution  of  aggregates  and  reasonable  snowfall 
rates . 

4.2.4  Aggregate  Riming 

The  growth  of  a  single  aggregated  particle  by  riming  can  be 
approximated  as 


dM  1 


E(a/c)  p  r 
o  c 


(4.17) 


Assuming  V  )>  V  and  applying  (4.17) 
a  c 

aggregates. 


to  a  population  of 


°°  dM  ' 

RM  =  —  /  — N(D)  dD 
C“  PoodtjRM 


which  gives  us  after  integration  and  substitution  of  (4.17)  and  (4.5) 


RM 

ca 


2.42 


1/2 


E( a/c ) 


r 

a 


r 

c 


(4.19) 


where  E(a/c)  is  evaluated  similar  to  individual  crystals  collecting 
cloud  droplets  (see  Cotton  et  al ■ .  1982). 
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4.2.5  Aggregate  Collision  with  Raindrops 
Aggregates  can  also  collide  with  supercooled  raindrops  leading  to  a 
loss  in  concentration  of  raindrops 

CD 

dglRjdR  =  _  y  n(R  +  d/2)2  |v(r)  _  v  (D) |E(R/a)^(R)N(D)dRdD 
dt  a 

0  (4.20) 


as  well  as  a  loss  in  concentration  of  aggregates 


dD  =  -  /  n(R  +  D/2)2|V(R)-V  (D)  |  E(R/a )  f»  (R)N(D)  dRdD 


and  providing  a  source  of  granpel.  The  change  in  mixing  ratio  of 
graupel  due  to  collision  between  aggregates  and  supercooled  raindrops  is 
given  by 


CL 


rag 


cl  -  J-  7 


arg 


o  o 


dt 


-dR  +  —  /  m(P)^-Pj-dD 

p  dt 

*o  o 


(4.22) 


which  after  some  manipulation  become 


CL  =  2. 53 Op  P 
rag  o> 


>R  ~  vjE(R/a> 


r  N  R3D  '24 
a  R  m  m 


[40R2+8R  D  +D21 
L  m  m  m  mJ 


(4.23) 


and 


a.  =6.32  r  Njv_-V  |E(R/a)(RZ  +  1.7  R  D 

arg  a  R'  R  a  a  m  n 


+  1.87  D  ) 

m 


(4.24) 


where  the  appropriate  raindrop  parameters  are  given  in  Cotton  et  al ■ 
(1982)  and  Tripoli  and  Cotton  (1980). 


The  rate  of  mass  growth  of  a  single  aggregate  particle  by  vapor 


deposition  can  be  approximated  as: 


=  4 »tC  G(T,D)  f (R  ) 
e 


(S.-l) 

1 


(4.25) 


where  C,  the  capacity  is  evaluated  as  a  spherical  particle  thus 
C  =  D/2, 


f(R  )  is  the  ventilation  function  and  S.  is  the  saturation  ratio  with 
e  i 

respect  to  ice.  The  change  in  aggregate  mass  due  to  vapor  deposition  is 
then  given  as 


®  .  dM  (D) 

VD  =  /  — - —  N(D)  dD 

va  p  dt 

A  r  A 


(4.26) 


which  when  integrated  over  the  entire  spectrum  of  aggregates  becomes: 

VD  =  4.03  D~1-4  p"1  G(T,P)  f(R  )  (S.-l)  r  (4.27) 

va  m  l  e  i  a 

4.2.7  Melting  of  Aggregates 

We  treat  the  melting  of  aggregates  in  a  similar  manner  to  the 
melting  of  graupel  particles  (see  Cotton  et  al . .  1982).  Consider  an 
aggregate  falling  through  a  cloudy  environment  at  temperature  T  >  0°C 
(T^)  and  accreting  cloud  droplets  having  a  temperature  T  and  colliding 
with  raindrops  also  at  temperature  T.  It  is  assumed  that  the  melting 
aggregate  is  in  thermodynamic  equilibrium  with  a  surface  temperature  T^ 
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dm 

L  — ® 

li  dt 


=  -2nD[Kf(R  )  (T-T  )  +  L  ^  f  (R  )  (p  -p  (TJ)1 

,  (  e  f  vi  e  rv  rvs  f  | 

5 1  t  L  J 

+  ^a|  ]r  ~  x  <4.28) 

I  dt  JCL  J 


dm 

_ ! 

dt 


c  (T-T.) 

w  f 


Application  of  (4.28)  to  the  spectrum  of  aggregates  give  us: 


ML  =  -4.03  f (R  )  L  1,4r  (K(T-T  )  +  L  .  <Pp  (r  -r  (T,)) 

ar  e  ilrl  m  a  f  vi  ro  v  vs  f 


-  r1-  (RM  -  CL  )  c  (T-T  ) 
L^  ca  rag  w  f 


(4.29) 


4.2.8  Aggregate  Collection  bv  Graunel 

Aggregates  can  be  collected  by  graupel  particles  thereby 

contributing  to  the  mixing  ratio  of  graupel  particles.  The  change  in 

concentration  of  aggregates  of  diameter  D  +  6D  by  colliding  with  graupel 

of  diameter  D  +  6D  is  given  by 
8  8 

.  .  «  n(D+D  )2 

SMJLLdD  =  -  /  - | V(D)-V(D  )  |E(a/g)N(D)N(D  )dDdD 

dt  „  4  '  8  8  8 

0  (4.30) 


The  change  in  mixing  ratio  of  aggregates  is  then 


CL  =  -  —  /  M(D) 


ag 


dN(D) 

dt 


dD 


(4.31) 


which  after  integrating  over  the  graupel  and  aggregate  spectra  becomes: 

CL  =  -1.57|\T-V-|E(a/g)N  r  (7.5  D2  +  3 .4A_1  D  +  A~2) 

ag  a  g*  6  tg  a  m  m 


(4.32) 
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4.2.9  Conversion  of  Aggregates  to  Graupel  Doe  to  Heavy  Riminp 
If  the  riming  rate  on  aggregates  exceeds  a  certain  amount  we 
presume  that  the  aggregates  become  a  source  for  graupel. 

Thus,  we  convert  aggregates  to  graupel  at  the  rate: 


CN  =  -  MAX  [RM  -  B] 
ag  ca 


(4.33) 


where  B  is  a  threshold  for  conversion.  We  estimate  B  as  20%  of  the 
graupel-equivalent  riming  rate  given  by 


RM  =  K  E(glc)  r  r5/6 
rg  c  a 


where  K  =  1.16  (g/C  ) 
rg  D 


M  N1 

W  1 


4.2.10  Summary  of  Model  Equations  for  Water  Snbstance 
Ihe  following  is  a  summary  of  the  relevant  equations  pertaining  to 
the  aggregate  population: 


—  fPl  8 

V  =  2.49  r 

8  Lpo  S 


(4.34) 


CL.  =  0.503  E(a/i)D_2'4p  |v  -V3|(2DZ+2D  D.  +  DZ)r  r.  (4.35) 

ia  m  o'  a  i  mmi  lai 


RM  =  2 
ca 
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E(a/c)  r  r 
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(4.36) 


(40  R2  +  8  R  D  +  D2)  _  _ 

- »  "■ - E(»/.>  po  |Vt  -  VRlriIr 

P1  m 


(4.37) 
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CL  =  0.25 
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CN  =  -MAX(Rm  -  B,0) 
ag  ca 

When  combined  with  the  water  substance  equations  described  in 
Cotton  et  al .  (1982)  a  complete  set  of  equations  governing  the 
distribution  of  water  substance  in  the  model  is  then  given. 

4.3  Summary  of  Conservation  Equations  for  Ice  Crystals 


As  can  be  seen  in  Eq.  (4.1),  the  rate  of  production  of  aggregates 


8  3 


is  proportional  to  N^.  Therefore  it  is  necessary  to  develop  a  parallel 

equation  for  the  continuity  of  ice  crystal  concentration. 

The  conservation  equations  for  concentration  of  ice  crystals  (N^) 

contains  the  following  processes: 

i)  nucleation  of  ice  by  deposition 

ii)  nucleation  of  ice  by  contact  freezing 

iii)  ice  multiplication  by  the  Hall et t-Mossop  mechanism 

iv)  removal  of  N.  by 
1 

a)  melting 

b)  conversion  to  graupel 

c)  collection  by  rain,  graupel  or  aggregates 

d)  precipitation 

v)  transport  and  diffusion 

In  order  to  simplify  the  transport  and  diffusion  of  crystal 

concentration,  it  is  convenient  to  predict  specific  concentration 

(N./p  ).  The  conservation  equation  is  then 
1  o 


d(N./p  ) 
i  o 

at 


-  u. 

j 


3(N1/p  )  3U.(N./p  ) 

_ L-.0-  _ i _ i— 


dx . 


3x  , 


+  NPR . 


N. 

+  — —[ML.  +  (CN ,  -  CL  .)  +  CL.  1 

PQri[  iv  ig  ri  lgj 


+  NNUA  .  +  max(0,NNUB  ,  +  NNUC  .  +  NNUD  .) 
vi  vi  vi  vi 

+  NSP  . 
vi 


(4.43) 


where  NPR^  is  the  precipitation  tendency  of  specific  crystals 

concentration,  ML.  is  the  loss  of  crystal  mixing  ratio  to  melting,  CN 

iv  ig 

-  CL^  is  the  loss  of  crystal  mixing  ratio  to  conversion  to  graupel  and 


CL.^  is  the  loss  of  crystal  mixing  ratio  to  collection  by  granpel.  We 
consider  nucleation  by  deposition  (NNUA^)  and  contact  nncleation  due  to 
Brownian  diffusion  (NNUB^  ) ,  thermophoresis  (NNUC^)  and 
dif fusiophoresis  (NNUD^).  Finally,  production  of  specific 
concentration  by  Hal lett-Mossop  splintering  (NSFP^)  is  considered.  For 
a  description  of  the  terminology  see  Cotton  et  al .  (1982).  The 
precipitation  term  NPR.  is  given  by 


NPR .  = 
1 


dz 


(4.44) 


where  V  is  the  crystal  terminal  velocity.  The  terms  ML.  ,  CN.  ,  CL  . 

iv  lg  ri 

and  CL.  are  described  bv  Cotton  et  al.  (1982).  Next  we  will  describe 
rg  J  - 

the  nncleation  and  splintering  models. 

4.3.1  Nucleation  bv  Deposition 
Empirical  results  of  Fletcher  (1962)  show  that,  for  a  water 
saturated  parcel,  depositional  nucleation  will  produce  a  relationship 
between  temperature  and  concentration  given  by: 


N.  =  N. 

1  10 


(4.45) 


where  N.  =10  *cm  ^ 
io 

Differentiating  Eq. 


and  p,  =  0.6C  and  T  =  (273.16-T), 
•  s 

(4.45)  with  respect  to  z  we  obtain 


3N.  0,T  3T 

TT  =  P2  Nio  *  Tz 


(4.46) 


If  we  neglct  the  variation  of  temperature  in  the  horizontal  and  time 


compared  to  the  vertical,  then 
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3N.1 

4,P 


WP2N 


io 


OT 
_ s 

d  z 


or,  for  specific  concentration 


NNUA  . 
vi 


max 


wp  N.  p.T 
Fx  10  s 

0.  , -  e 


3z 


(4.47) 


(4.48) 


4.3.2  Contact  Nucleation 

For  a  model  of  contact  nucleation  we  shall  follow  the  model  that 
Young  (1974)  referred  to  as  model  A.  Young  considers  contact  by 
Brownian  diffusion,  thermophoresis  and  dif fusiophoresis .  Brownian- 
diftusion  contact  nucleation  results  from  the  random  collision  of 
aerosol  particles  with  cloud  droplets.  Thermophoresis  contact 
nucleation  occurs  due  to  the  attraction  or  repulsion  of  aerosol 
particles  to  the  droplet  along  the  gradient  of  thermal  diffusion. 

Dif fusiophoresis  contact  nucleation,  on  the  other  hand,  occurs  due  to 
the  attraction  or  repulsion  of  aerosol  particles  to  the  droplet  along 
the  gradient  of  vapor  diffusion.  Hence,  the  net  effect  of 
thermophoresis  and  dif fusiophoresis  will  be  to  inhibit  contact 
nucleation  in  regions  of  supersaturation  and  have  the  opposite  effect 
where  subsaturation  occurs.  Since,  as  Young  points  out,  thermophoresis 
is  dominant,  these  processes  can  be  important  to  the  problem  of  aircraft 
icing . 

For  Brownian  diffusion,  the  change  in  specific  concentration  is 

NNUB  =  F,  DF  (4.49) 

ri  1  ar 


where  DF  is  the  aerosol  diffusivity  and  where 

ar 
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4nR  N  N 
c  c  t 


(4.50) 


We  approximate  the  activated  aerosol  concentration  by 


N  =  N  (270.16  -  T  ) ' 
a  ao  c 


(4.51) 


where  N  =  2  x  10  1  cm  3  at  sea  level.  Young  assumed  N  decreased 
ao  ao 

linearly  with  height  to  1  x  10  ^  cm  3  at  5000m  MSL.  We  will  assume  a 

-1  -3 

constant  value  of  2  x  10  cm  for  now.  Given  the  cloud  droplet 

concentration  (N  ),  and  the  cloud  water  mixing  ration  r  .  the  cloud 
c  c 

droplet  radius  is  given  by 


3nNcpwJ 


(4.52) 


Similarly,  for  thermophore  tic  contact  nucleation,  we  scale  Young's 
equations  to  obtain: 


NNUCvi  -  F,  F2  ft 


(4.53) 


where 


F  =  -  (T-T  ) 
2  P  C 


(4.54) 


and  K  is  the  thermal  conductivity  of  the  air,  T  is  the  cloud  droplet 

c 

temperature  and  p  is  the  atmospheric  pressure.  The  function  f^  is  given 


(4.55) 


•4[1+1.45E  +  .4E  exp (-1/E  )](E+2.5E  E  ) 

_ n _ e _ l _ a _ a.  a 

(1+3E  ) (2E+5E  E  +E  ) 
n  ana 


where  E  is  the  Enudsen  number  given  by 
n 


_  LJ7.T 

n  2  88pr  ^ 


(4.56) 


and  r  is  the  assnmed  aerosol  radios  of  3  x  10  cm.  The  aerosol 

a 

-4  -1  -1  -1 

thermal  conductivity,  E^,  is  assumed  to  be  5.39  x  10  ergs  cm  E 
Finally,  we  approximate  diffusiophoresis  contact  nucleation  by 


(e  -  e  (r  ))  (4.57) 

NNUD  =  F„  DF  — J - S — 

v  i  1  v  p 


where  DF^  is  the  vapor  diffosivity  and  e^  and  e  (r^)  are  the  vapor 
pressures  at  infinity  and  the  droplet  surface.  For  a  droplet  in  thermal 
equilibrium  it  can  be  shown  that: 


DF 

v 


3  -  e(r  )) 

_v _ g — 

P 


R  T 


(4.58) 


therefore  we  can  express  diffusiophoresis  contact  nucleation  as 


R  T 


NNUDbi  "  F1F2  L 


(4.59) 


vl 


Bence  the  sum  of  the  contact  nucleation  may  be  written 


NNUB  .  +  NNUC  .  +  NNUD  ,  =  F,  DF  +F-(f 
vi  vi  vi  1  ar  2 


F  RT  l 

*  DF  +F-(f  -  r*-) 

l  11  2  L1V  J 


(4.60) 


Since  the  quantity  F  is  dependent  on  T-T  ,  the  proper  evaluation 

z  c 


A 
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of  thermophoretic  and  diffusiophoretic  contact  nacleation  must  include  a 
realiatic  eatimation  of  thia  temperature  difference.  Assuming  thermal 
equilibrium, 

G(T,p) (S-l)L  (4.61) 

T-T  =  -  - 

c  K 


where  G(T,p)  is  the  portion  of  mass  tendency  of  a  water  particle 
resulting  from  vapor  diffusion  to  or  from  the  particle,  since  we 
normally  diagnose  cloud  water  in  the  model  (S-l)  is  zero  thus  both 
thermophoresis  and  dif fusionphoresis  would  always  be  zero.  In  order  to 
estimate  a  supersaturation,  we  assume  (S-l)  results  from  a  balance 
between  production  and  removal  mechanisms  of  supersaturation.  We  assume 
the  dominant  removal  is  by  diffusional  growth  of  cloud  droplets  and  ice 
crystals  and  neglect  all  other  growth  processes.  The  dominant 
production  is  by  temperature  change  in  response  to  vertical  motion. 
Following  Tripoli  and  Cotton  (1980)  we  write  the  production  as 


P 

s 


wr  eL  , 

—5  *  (—Si  _j) 

Tv  R  CT 
_ E _ 


eL2  r 

_ yj . s 

RC  T2 
P 


+  1 


(4.62) 


where  w  is  vertical  velocity,  rs  is  saturation  mixing  ratio,  Ty  is 

virtual  temperature,  g  is  gravity  acceleration,  8  =  .611,  and  C  is 

P 

specific  heat  at  constant  pressure.  The  dissipation  of  supersaturation 
Dg  is  given  by 
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r 

r 

v 


dr  dr . 

— &  + 
dt  dt 


where 


vd 


—  G(T,p) (S-l )N  R 
Po  c  c 


(4.63) 


and 


=  VD 


vi 


(4.64) 


For  a  steady  state  anpersataration. 


P 

s 


(4.65) 


r 

or  coartiining  Eqs.  (4.63)  -  (4.65)  and  (4.52)  and  assuaiing  (~i)  ~  1 

v 

obtain 


we 


(S-l) 


-p  (VD  J 
_S _ Xi_ 


dr 
_ 8 

AL 


4nG(T,p)N  R 
c  c 


(4.66) 


4.3.3  Ice  Multiplication 

To  aiodel  ice  multiplication  we  shall  make  use  of  a  parameterization 
of  the  Hallett-Mossop  theory  developed  by  Gordon  and  Marwitz  (1981). 

Two  mechanisms  are  parameterized: 

1.  Hallett  and  Mossop  (1976)  reported  that  approximately  350  ice 

_3 

splinters  are  produced  for  every  10  grams  of  rime  accreted  on 
graupel  particles  at  -5°C.  The  parameterization  evaluates  the  rate 
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of  crystal  production  as 


P  (I)  =  3.5  x  It)5*'1  J7I  .  .  f  (T  ) 
p  dtJ riming  1  f 


(4.67) 


dMT 

where  dtJriming  *s  t*ie  rtoiag  rate  of  an  ice  crystal  or  grapel 
particle.  The  function  f(T^)  is  given  by: 


IT  -  268.16 


W  IT  -  268.16 


T  >  270.16 


270.16  >  T  >  268.16 


268.16  >  T  >  265.16 


(4.68) 


265.16  >  T 


where  T^  is  the  temperature  of  the  crystal. 

2.  Approximately  one  ice  crystal  is  produced  for  every  250  drops 

-4 

greater  than  2  x  10  cm  radius  accreted  onto  each  graupel  particle 
at  -5°C  (Nossop,  1976) .  The  ice  crystal  production  rate  per 
particle  is  then 


V11’  ■  isS  -?  V'^u’Vi'V 


(4.69) 


where  is  the  concentration  of  cloud  droplets  greater  than  12  x 
-4 

10  cm  in  radius,  Efg/c^j)  is  the  efficiency  of  collision  between 

-4  — 

graupel  particles  and  cloud  droplets  greater  than  12  x  10  cm,  V 

8 

is  the  terminal  velocity  of  graupel  and  D  is  the  diameter  of 

8 

graupel . 

To  estimate  the  concentration  of  droplets  greater  than  12  x  10  cm 
in  radius,  we  assume  that  the  cloud  drovlet  distribution  is  given  by  a 
Gamma  distribution  of  the  form: 


AD- A 142  690 
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f  ( r) 


M  a-1  -r/p 
Nr  e  r 
c _ 

n  a)  pa 


0  <  r  <  00 


(4.70) 


where  Cotton  (1970)  showed  that 


a 


(4.71) 


where  is  the  radius  dispersion.  We  choose  the  typical  value  of  y = 
0.18. 


Furthermore,  Cotton  showed  that 


0  = 


3M 


[4na(a+l ) (a+2) 


1/2 


(4.72) 


where 


M 
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r  P 
co 


N 

c 


(4.73) 
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The  concentration  of  droplets  having  radii  greater  than  12  x  10  cm  is 
then 


Nj2  =  /  f (r)dr 

r12 


c  r  a-1  -r/p 
-  J  r  e  dr 

r(a)p*  o 


(4.74) 


Letting  X  =  r/0(M  ),  we  find: 


N.,  =  N  f,(M  ) 
12  c  3  c 


(4.75) 


/ 

A 


where 
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f,(M  ) 
3  c 
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~r12/KMc) 


va-l  -X 
X  e  dX 


(4.76) 


The  evaluation  of  this  incomplete  gamma  function  is  shown  in  Figure  1. 
It  is  also  shown  to  a  good  approximation  that 


f,  (m  ) 
3  c 


0  m  <  1.26  x  10'6g 

c 

+2.27 In  m  +  13.39  1  .26  x  10_6g<  m  <  3.55  x  10_6g 

C  c  — 

1  a  >  3.55  i  10 

c 


(4.77) 


Applying  Eqs .  (4.69)  and  (4.67)  to  simulated  riming  tendencies  we  obtain 


NSP 


ri  pr 


N 

p.p  +  4x10  3  irf-O*  ) 
5ro  3  c 
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ff,(T.)CL  .  +  f, (T  )CL  +  f, (T  )CL  1  (4.78) 

[li  ci  1  g  eg  la  caj 


where  =  3.5  x  10^g  * 


Eq.  (4.78)  is  then  complete  except  for  the  evaluation  of  T. ,  T  and  T  . 

i  g  a 


Assuming  thermal  equilibrium. 
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T  =  T  +  a*** 
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Approximation  of  incomplete  gamma  function  to  determine 
percent  of  clond  droplets  size  greater  than  10  pm  assum 
droplets  in  gamma  distribution,  dispersion:  y  =  .018. 
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T  =  T  +  0 .25p  —  (L  VD  -  L..CL  ) 

•  Kw  —  vi  va  il  ca 

a 


(4.81) 


where  is  the  mean  ice  crystal  diameter  and  N  is  the  graupel 
concentration. 

4.4  Results  of  Experiments 

The  application  of  our  aggregation  and  ice  crystal  concentration 
models  to  an  orographic  case  study  has  thus  far  demonstrated  the 
fol lowing 

1)  The  ice  nucleation  deposition  equation  produced  too  many  crystals  at 
cold  temperatures  and  was  insignificant  compared  to  other  terms  at 
warm  temperatures.  Therefore  we  elected  to  ignore  it  at  the 
present . 

2)  Aggregates  were  produced  only  when  the  collection  efficiency  was 
enhanced  at  the  colder  temperature,  i.e.,  -12  to  -15°C. 

3)  Supercooled  water  diminished  as  nucleation  took  place.  The  mixing 
ratio  of  supercooled  water  was  found  to  be  very  sensitive  to  contact 
nucleation.  Apparently  a  close  balance  occurs  between 
thermophoresis  and  diffusiophoresis  contact  nucleation.  When  one  is 
limited,  many  more  crystals  are  produced  and  much  less  liquid  water 
exists.  Th i s  is  because  a  larger  number  of  small  crystals  drive  the 
Bergeron-Findeisen  growth  process  more  efficiently. 

More  experimentation  is  being  done,  but  so  far  we  have  found  an 
interesting  relationship  between  nucleation  and  the  existence  of 
supercooled  liquid  water.  The  production  of  aggregates  also  affected 
the  existence  of  liquid  water  because  it  removed  large  numbers  of  ice 
crystals  and  thus  affected  the  overall  balance. 
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5.0  Inclusion  of  Ice  Processes  in  the  Existing  ID  Cloud/Turbulence 
Model 

5.1  Introduction 

Chen  and  Cotton  (1983b)  have  described  a  one-dimensional  higher 

order  turbulence  model  which  was  applied  to  the  simulation  of  a  cloud- 

capped  atmospheric  boundary  layer  (ABL) .  The  model  has  exhibited 

considerable  skill  in  quantitatively  reproducing  the  features  of  a 

marine  ABL  that  have  been  predicted  with  Deardorff’s  (1980)  three- 

dimensional  large-eddy  simulation  (LES)  model.  More  recently,  the  model 

has  been  tested  against  observations  of  a  stable,  marine  ABL  reported  by 

Brost  et  al .  (1982ao) .  The  model  has  also  been  successful  in 

reproducing  the  observed  features  reported  by  Brost  et  al .  (1982ab) 

including  the  effects  of  wind  shear  and  of  drizzle  on  cloud 

organization.  It,  therefore,  seemed  natural  to  investigate  the 

feasibility  of  expanding  this  cloud  system  into  the  ice-phase.  To  begin 

with  we  shall  adopt  the  following  notation. 

Mean  Equations 

VD  .  s  vapor  diffusion 
ab 

CL  .  -  collection 

ab 

CN  ,  =  conversion 

ab 

ML  ,  s  melting 
ab 

NUA  =  nucleation  by  deposition 

ab 

NOB  .  =  nucleation  by  Brownian  collection 

ab 

NBC  .  =  nucleation  by  thermophoresis 

ab 

NUD  =  nucleation  by  dif fusiophores is 

ab 

SP  =  splintering  formation  of  ice  crystals 

ab 


=  shedding 


Subscripts  of  the  above  terns  define  the  source  (second  subscript)  and 
sink  (first  subscript)  water  categories  between  which  the  transfer  is 
made.  The  subscripts  are  also  used  with  mixing  ratios  to  define  which 
category  the  mixing  ratio  is  assigned.  The  water  categories  are  as 
fol lows : 

v  s  vapor 
c  ®  cloud  water 
r  —  rain 
i  a  ice  crystals 
g  £  graupel 
a  s  aggregates 

When  a  process  description  is  preceeded  by  the  letter  'N'  the  term 
refers  to  an  ice  crystal  concentration  tendency. 


Turbulence  Equations 


VD'  ' 

ab 

5 

turbulence 

contribution 

to 

vapor  diffusion 

CL' ' 

ab 

s 

turbulence 

contribution 

to 

col lection 

CN ' ' 

ab 

wm 

turbulence 

contribution 

to 

auto  conversion 

ML’ ' 
ab 

m 

turbulence 

contribution 

to 

melting 

SH' ' 

ab 

= 

turbulence 

contribution 

to 

shedding 

The  turbulence  contributions  to  the  nucleation  by  deposition,  Brownian 
collection,  thermophoresis,  dif fusiophoresis  and  splintering  of  ice 
crystals  are  not  included. 

rT  -  total  water  mixing  ratio 

0^  =  ice-liquid  water  potential  temperature 

The  method  of  deriving  ensemble-average  mean  equations  and  second  order 
turbulence  equations  has  been  reported  by  Cotton  and  Anthes  (1983). 
Basically,  we  decompose  every  variable  p  into  Reynold’s  averaging  mean  p 
and  perturbation  p' ' ,  i.e.  p  =  p  +  f' ' *  In  the  following  sections,  we 
will  introduce  the  summary  of  the  mean  equations,  then  the  perturbation 
equations . 

5.2  Summary  of  the  Mean  Equations  and  Micronhvsical  Processes 


5.2.1  Mean  Equations 

Rain  Water 
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5.2.2  Microphysics 

The  parameterization  of  microphysics  follows  the  same  scheme 

described  in  Section  4.  However,  because  of  the  introduction  of 

turbulence,  some  additional  terms  can  be  found  in  the  following 

equations.  For  instance,  the  rate  of  accretion  of  cloud  water  by  rain 

water  is  given  by  CL  =  a_ r  r  (Note:  Coefficients  a,,  a«  ...  are 

cr  2  c  r  12 

defined  in  the  Appendix,  Section  5.5).  Here  the  average  accretion  rate 

is  given  by  CL  =  a„  (r  r  +  r  *'r  '*).  Chen  and  Cotton  (1983)  found 
cr  2  c  r  c  r 

that  the  correlation  term  r  '*r  ’*  cannot  be  neglected  in  a 

c  r 

s t ra tocumul us  environment.  Both  terms  r  r  and  r  ' ’r  ' '  have  the  same 

c  r  c  r 

order  of  magnitude. 

The  parameterization  of  nucleation  by  deposition,  nucleation  by 
Brownian  collection,  nucleation  by  thermophoresi s ,  nucleation  by 
dif fusiophoresis  and  splintering  of  ice  crystals  are  the  same  as 
described  in  Section  4.  The  turbulence  contributions  to  the  above 
mentioned  processes  are  not  considered  here. 
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Collection 

The  parameterization  for  the  mean  accretion  rate  of  rain,  granpel, 
aggregates,  ice  water  and  ice  crystal  can  be  summarized  as: 


‘25 


(r  r  +  t  '  ' t  ' ' ) 

a  g  a  g 


(5.6) 


CL 


ca 


=  a,,  (r  r 
31  c  a 


+  r 


(5.7) 
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(r  r  +  r  "r  ") 
eg  eg 
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CL 


ci 
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(N.r  +  N.'  'r.")  h(T  ) 
i  c  11  s 


(5.9) 


CL  =  a„  (r  r  +  r  ’  'r  '  ') 
cr  2  c  r  c  r 


(5.10) 


CL  =  a,,  (r  r 
ra  11  r  a 


+  r 


(5.11) 


CL 
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(rirr 


+  r. 


(5.12) 
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(r  r  +  r  '  *r  '  *) 

r  *  r  8 


(5.13) 
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CL. 
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"r.”)  h(T  ) 
l  s 


(5.14) 


CL. 
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20 


(r.r  +  r  ."r  ")  h(T  ) 

i  g  x  g  s 


(5.15) 


Conversion 


The  mean  conversion  can  be  written  as: 


CN 


ag 


=  (i 


_^(r  r 
26  c  a 


+  r 


)  +  a27  <rari 


r  ”))  h(T  ) 
r  s 


(5.16) 


CN  =  a,(r  r  h(r  -  r  ) 

cr  3c  c  c  cm 


(5.17) 


CN. 


=  a  , (N. r . 

22  li 


+  N. 


') 


(5.18) 


N. 

x. 


CN  =  (max  ( a,  fi(  r  r.  +  r  *’r.")  -  — C  ,0)+a«o(r  r.+r  "r."))k(T  ) 
ig  28  ci  ci  pm  29riri  s 


(5.19) 


Melt ina 


The  melting  rate  can  be  represented  by 


10} 


ML  =  U,,r  +  a,,(r  r  +  r  "r  **)  +  a,.(r  r  +  r  "r  "))h(-T  ) 
ar  12  a  13  a  c  a  c  14  a  c  a  c  s 


(5.20) 


ML  =  (a,r  +  a_(r  r  +  r  "t  ")  +  a_(r  r  +  r  "r  "))h(-T  ) 
gr  6  g  7  eg  eg  8  r  g  r  g  s 


(5.21) 


ML.  =  a. _  r.  h(-T  ) 
iv  18  l  s 


(5.22) 


Nncleation 

The  mean  nncleation  rate  can  be  summarized  as: 


NCL .  =  a  (N.  r  +  N.'  'r  '  ') 

la  38  l  a  i  a 


(5.23) 


NCL.  =  a,,(r  N.  +  r  "N.”) 
ig  36  g  i  g  l 


(5.24) 


NCNi.  ■  *37<Ni2  *  "l  > 


(5.25) 


NCN.  =  (a,,(r  N.  +  r  "N.")  +  a.,(r  N.  +  r  ”N."))h(-T  ) 
ig  34  c  l  c  i  35  r  i  r  i  s 


(5.26) 
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3T 

exp ( 8 -T  )  w  # v •  *»\*  r 

2  s  d  z  J  s 


L,0.1  h(T  ) 


(5.27) 
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_  4nR  N  n  m 
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(5.30) 


NML .  =  a.  N. 

iv  33  l 


(5.31) 


Shedding 

The  mean  shedding  terms  are: 


x 

A 


(a,,(r  r  +r  "r  ")  +  a,,(r  r  +r  "r  "))h(-T  ) 
1 5  a  c  a  c  16  a  r  a  r  s 


(5.32) 


SH  =  ( a.(r  r  +r  "r  ")  +  a1A(r  r  +r  "r  "))h(-T  ) 
gr  9  c  g  c  g  lOrgrg  s 


(5.33) 


The  mean  splintering  term  is: 


3  N 

„-3  c 


SP  .  =  ~  p.p  +  4  i  10  r  f, (m  )  [f, (T.)CL  .  +  f. (T  )CL  +  f. (T  )CL  1 

vi  p  5  o  3c  |  1  l  ci  1  g  eg  la  caj 

c 


(5.34) 


where  CL  . ,  CL  and  CL  are  mean  accretion  terms  which  can  be  found  in 
ci  eg  ca 

the  section  on  collection. 


The  mean  vapor  deposition  terms  are: 


VD  =  a,A  r 
va  30  a 


(5.35) 


VD  =  a„„  r 
vg  23  g 


(5.361 


^vi  =  *17  Ni  h(V 


(5.37) 


106 


VD 

vr 


al 


r 

r 


(5.38) 


5 .3  Summary  of  Turbolence  Flaxes  and  Covariances 
5.3.1  Introduction 

The  procedure  to  derive  the  formulation  for  the  mean  variables, 
turbulence  fluxes  and  its  covariances  follows  the  same  methodology  as 
described  by  Manton  and  Cotton  (1977).  Since  the  introduction  of  rain, 
ice,  graupel,  aggregate  and  ice  crystal  in  this  report,  the  equations 
for  the  turbulence  fluxes  and  covariances  are  much  more  complicated  than 
the  formulation  described  by  Chen  and  Cotton  (1982,  1983ab).  In  spite 
of  the  complication  of  the  equations,  some  simple  and  consistent 
patterns  of  equations  can  be  derived.  In  the  next  section  the 
terminology  for  the  correlation  between  various  terms  is  introduced. 

The  summary  of  the  turbulent  microphysics  is  described  in  a  later 
section. 


5.3.2  Terminology 

In  the  equations  below  many  terms  have  similar  structure,  but  their 
difference  depends  upon  what  variables  are  taken  into  account.  The 
perturbation  variables  are  represented  by 


#'  = 


u-.  rT".  e:;. 


,  r  ",  r  ",  N."/p 

g  a  i  c 


The  variable  r  '  '  denotes  the  perturbation  mixing  ratios 


( r 


•  f  —  *  # 

#  1  • 


N.'  ’ 

r  "  and  ~i~) 

a  p 


The  following  equations  represent  the  production  terms  for 
turbulent  fluxes  (Equation  5.39)  and  turbulent  covariances  (Equation 
5.40).  The  various  terms  on  the  right  hand  side  (RHS)  of  Equation 
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(5.39)  represent:  shear  prodnction,  advection,  eddy  transport,  isotropy 
redistribution  and  buoyancy.  The  terms  on  the  RHS  of  Equation  (5.40) 
represent:  advection,  shear  productions,  eddy  transport  and  dissipation 

of  the  covariances. 


A(r  ")  = 
T 


-  u  ’] 
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a  - 
a*,  “i 
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-  a  — 

U  .  -  u.* 

j  3x  .  i 
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dx  . 
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(5.39) 
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a 

dx  . 
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u.'  *#’  'r  '  ' 
J  r 


-2  e(f  1  ,r  ") 
r 


(5.40) 


where  e  (f'',  r^'  ’)  is  the  dissipation  term  which  can  be  parameterized 

as 


e(*",r  ’  ')  = 
r 


M#  * '  rr*  ’ 


(5.41) 


where  p  =  0.31,  rA  is  defined  in  Chen  and  Cotton  (1983a). 

In  the  following  equations,  the  turbulent  microphysical  processes 

to  produce  rain,  aggregate,  graupel,  ice  and  ice  crystal  are  denoted  by 

R’’,  A*’,  G  '  ' ,  I'*,  and  IC  ' ' .  The  mean  equations  for  those  processes 
p  p  P  p  P 

correlated  with  other  turbulence  variable  (f’)  can  be  written  as 
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(5.64) 


(5.65) 


(5.66) 


(5.67) 


(5.68) 


(5.69) 


(5.70) 


(5.71) 


(5.72) 
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Turbulent  flox  of  ice  crystals 
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5.3.4  Definitions  of  Turbulent  Microphvsics  Terns 
In  this  section,  the  terns  on  the  right  hand  side  of  Equations 
(5.43),  (5.44),  (5.45),  (5.46),  and  (5.47)  are  defined  according  to  the 
following  processes:  collection,  conversion,  melting,  shedding,  and 
deposition. 

Collection 
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Melting 
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Splintering 
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Not  parameterized 
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Nucleation  by  deposition.  Brownian  collection 
thermophoresis,  dif fusiophoresis 

Since  very  little  is  known,  especially  for  the  turbulence 
contribution  to  the  nucleation  process  by  deposition.  Brownian 
collection,  thermophoresis  and  dif fusiophoresis  is  the  incorporation  of 
turbulent  fluctuations  in  the  parameterization  processes  would  be 
meaningless.  Thus,  the  parameterization  for  the  turbulent  related 
nucleation  is  not  proposed  in  this  report. 


1  17 


APPENDIX 


Suaunary  of  coefficients 
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5 .4  Simplified  ID  Cloud/Turbulencc  Model 

As  one  can  readily  see,  expansion  of  the  higher-ordered  turbulence 
model  to  the  ice-phase  results  in  an  extremely  complex  system  of 
equations.  Application  of  this  system  to  a  one-dimensional  model  has 
possibilities,  but  the  extension  to  the  2D/3D  flow  model  appears 
hopeless.  We  have,  therefore,  decided  to  develop  a  simpler  "hybrid" 
version  of  Chen  and  Cotton's  (1983b)  cloud/turbulence  model.  Following 
Mellor  and  Yamada's  (1974)  terminology  we  will  refer  to  Chen  and 
Cotton's  model  as  a  level  3.5  turbulence  model  with  the  new  form  called 
a  level  2.5  model.  The  objective  is  to  develop  a  turbulence  model  in 
which  the  number  of  prognostic  variables  is  greatly  reduced,  yet 
essential  features  of  the  level  3.5  model,  such  as  cloud  fractional 
coverage  and  entrainment  rates  are  properly  predicted.  The  level  2.5 
model  will  be  tested  and  developed  against  the  level  3.5  model,  first  in 
dry  ABL  applications  and  then  a  liquid-phase-only  cloudy  ABL.  After  we 
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have  decided  we  have  a  satisfactory  "compromise’'  model,  it  will  then 

be  extended  to  the  ice-phase  system  described  in  the  previons  sections. 

The  proposed  level  2.S  model  has  three  prognostic  equations 

(q2.  a' ’2)  and  diagnostic  equations  for  the  remaining  turbulence 

variables.  The  diagnostic  equations  compute  uV,  uj'a’’,  and  a'  ’b*  '  . 

In  the  derivation  of  diagnostic  equations,  we  assume  no  time  derivative, 

no  vertical  divergence  of  the  third  order  term  and  the  dissipation  of 

the  turbulence  energy  (e)  is  balanced  by  shear  production  (P  )  and 

s 

buoyancy  production  (P  )  of  the  turbulence  energy.  The  summary  of 

b 

equations  for  the  level  2.5  model  is  given  as  follows: 

r~q2  =  -  /-(u"2w'  '  +  v-  ’V  '  +  V77  +  2  —  P"w") 

3t’  dz  p 


+  2  (P  +  P  -  e) 
s  b 


(5.106) 


where 


-  u'  'w*  '  ~  u  -  v*  'w* '  ~~  v 
dz  dz 


(5.107) 


pb  =  -  R1  V'9U  -  R2  *"*' 


(5.108) 


e  =  q  /t 


(5.109) 
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KH  d  z  rT 


(5.121) 


Substitution  of  Equations  (5.120)  and  (5,121)  into  Equation  (5.117)  and 
(5.108)  results  in 
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ilrT 


(5.122) 


Pb  =  CKH 


(5.123) 


where 
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2  dz  V 
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(5.126) 


Substitution  of  Equations  (5.118),  (5.119)  and  (5.123)  into  Equation 
(5.112)  results  in 


+  D  K„  +  E 


(5.126) 


where 


2t  , ,  d  ~.2  ,  d  2 . 

=  -  — 1 ((  —  u)  +  (r~  v)  ) 


3c,  dz 


dz 


(5.127) 


E  =  *  C 
3  C1 


(5.128) 


Substitution  of  Equations  (5.126), 


(5.121)  and  (5.122)  into  Equation 
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(5.116)  results  in 
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Equation  (5.130)  can  be  reduced  to 
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where 


A  =  /<a  e-,) 

il  3z  il 


(5.131) 


Substitution  of  Equations  (5.118),  (5.119),  (5.120)  and  (5.121)  into 
Equations  (5.114)  and  (5.115)  results  in 


=  f(km  +  v  ^«i 


(5.132) 
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where 
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Substitution  of  Equations  (5.126),  (5.132)  and  (5.133)  into  Equation 
(5.113)  results  in 


u!  'w* 
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■  ^'-<3  *  D  H  *  E  V  £  *1  -  V(*H*h)  f? i 


a  - 


a- 


-  *2  0<t»  *  V  fz  li 


(5.136) 


Equation  (5.136)  can  be  reduced  to 

*»  ‘  ^‘<3  I2  *  »  *K  *  E  V  *  V(tM  *  V 

*  *2  G  ,f«  *  V1 


(5.137) 


We  can  obtain  K  and  K  by  solving  Equations  (5.130)  and  (5.137). 

5 . 5  The  Parameterization  of  Dissipation  and  the  Turbulence  Time 
Scale 

As  mentioned  before,  e,  represent  the  mean  rates  of 

turbulence  energy  dissipation,  destruction  of  0 !  '  ^ ,  r  '  ,  and  9''r'' 

il  T  ilT 

respectively.  Various  ways  can  be  found  in  the  literature  to 
parameterize  those  dissipation  and  destruction  terms.  For  example, 

Zeman  and  Lumley  (1976)  use  prognostic  equations,  while  Andre  et  al . 
(1977)  and  Sun  and  Ogura  (1980)  try  to  diagnose  those  terms.  Both 
Andre  and  Sun  and  Ogura  approach  the  above  problems  by  adopting 
Blackadar's  (1962)  length  scale  which  is  very  popular  and  is  employed  by 
many  authors.  The  advantage  of  using  Blackadar's  length  scale  (shown  by 
Eq.  5.143)  is  its  simplicity.  Blackadar's  length  scale  tells  us  the 


mixing  length  near  the  snrface  is  proportional  to  kz ,  where  k  is  the 
Von-Karman's  constant.  This  mixing  length  can  be  derived  from  the 
similarity  theory.  As  in  the  well  mixed  layer,  Blackadar's  length  scale 
converges  to  a  constant  mixing  length  which  is  around  50  m  ~  60  m.  The 
disadvantage  in  using  Blackadar's  length  scale  is  that  it  may  not  be 
appropriate  to  apply  it  to  a  cloud  layer.  The  reason  can  be  explained 
as  follows:  The  release  of  latent  heat  due  to  the  condensation  or  the 
convective  instability  created  by  the  cloud  top  radiation  cooling  and 
the  cloud  base  radiation  warming  may  be  accompanied  by  the  production  of 
turbulence  kinetic  energy  (TKE).  Osusally,  the  dissipation  of  TKE 
should  respond  very  quickly  to  the  production  of  TKE.  Therefore, 
Blackadar's  length  scale  can  not  function  well  in  this  instance. 

Therefore,  a  new  scheme  to  diagnose  a  turbulence  length  scale  is 
proposed  in  this  study.  This  new  scheme  not  only  contains  Blackadar’s 
length  scale,  which  is  valid  in  the  unstable  ABL,  but  also  can  deal  with 
an  unstable  cloud  layer  as  well. 

This  new  turbulence  length  scale  is  defined  by 
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h>  =  *75 


(5.140) 
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Lg  =  kz / ( 1  +  kz/Az) 
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(5.142) 


(5.143) 

(5.144) 


Cj  =  .102 

The  length  scale  represents  the  turbulence  length  scale  for  the 
stable  stratified  layer.  The  Blackadar's  length  scale  is  denoted  by  Lg. 
The  buoyancy  production  of  turbulence  kinetic  energy  is  defined  by  Eq. 
(5.144).  We  define  Eq.  (5.142)  based  on  the  fact  that  the  buoyancy 
production  of  turbulence  kinetic  energy  (B)  is  always  a  dominant  term. 
Dissipation  therefore  should  always  adjust  to  B.  Eq.  (5.142), 
indicated  that  the  dissipation  is  about  60%  of  B.  Kaimal  et  al .  (1976) 
found  that  the  mid-layer  dissipation  rate  is  about  0.4  -0.5times  the 
buoyant  production  rate. 

Three  types  of  turbulence  time  scales  (x,  Xg,  r ^ )  have  been 
discussed.  These  are  defined  by 
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=  u  \fs  L  (5.145) 

clq 

Tq  =  0 . 5  t  (5.146) 

t3  =  t/3.  (5.147) 

Eq.  (5.146)  is  defined  by  assuming  that  t  /x  is  about  .5  as 
discussed  by  Zeman  and  Lumley  (1976). 

5 .6  Summary 

The  scheme  for  the  level  2.5  turbulence  model  is  summarized  in  this 
report.  We  hope  to  implement  this  scheme  into  the  ID  version  of  the  CSU 
cloud/mesoscale  model  in  the  near  future. 

Because  of  the  various  options  of  the  CSU  cloud/mesoscale  model,  we 
have  presented  an  example  to  show  the  relationship  between  the  ID  and  3D 
model  code.  The  newly  formulated  scheme  and  coefficients  can  be  tested 
in  the  ID  code  and  then  applied  to  the  3D  code.  This  interaction 
between  the  ID  code  and  the  3D  code  can  be  very  crucial  for  the  future 
development  of  the  3D  model. 
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6.0  Expansion  of  the  3D  Clond/Mesoscale  Model  into  a  Hydrostatic 

Mesoscale  Model 

A  hydrostatic  option  has  been  successfully  implemented  in  the  non¬ 
hydrostatic  3D  clond/mesoscale  model  (Tripoli  and  Cotton,  1982; 
hereafter  referred  to  as  TC).  It  provides  greater  computational 
efficiency  on  large  scale  simulations  (Ax  >  10  km).  Since  it  is  an 
option  to  the  cloud/mesosca le  model,  most  of  the  numerical  formulations 
are  the  same  as  described  in  TC.  The  differences  will  be  detailed  in 
the  following  sections. 

6 . 1  Features  of  the  Hydrostatic  Model 

Most  aspects  of  the  model  described  in  TC  are  the  same  including 
the  coordinate  system,  thermodynamics,  microphysics,  advective 
operators,  turbulence  parameterization  and  lateral  and  bottom  boundary 
conditions.  The  only  differences  arise  in  the  pressure  and  vertical 
velocity  computations.  Vertical  velocity  may  be  diagnosed  with  a  choice 
of  the  incompressible,  anelastic,  or  compressible  continuity  equations. 

A  compressible  or  incompressible  version  of  the  hydrostatic  option  is 
available  with  the  compressible  version  formulated  in  a  "time-split" 
mode.  This  allows  for  a  factor  of  2-3  increase  in  computational 
efficiency. 


6.1.1  Vertical  Velocity  Diagnosis 
The  vertical  velocity  may  be  diagnosed  with 
incompressible,  2)  anelastic,  or  3)  compressible 
(All  notation  after  TC.) 

1)  incompressible  continuity  equation 

9w  du  _  dv 

dz  dx  dy 


a  choice  of  the  1) 
continuity  equations. 


(6.1) 


j 
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2)  anelastic  continuity  equation 
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dp  w  dp  u  dp  v  (6.2) 

O  _  O  o 

dz  dx  dy 


3)  compressible  continuity  equation 

dp  w  dp  u  dp  v  ,  (6.3) 

ro  _  o  _  ro  dp 

dz  dx  dy  dt 


To  diagnose  the  compressible  term  p  is  substituted  from  the 
hydrostatic  equation  and  the  last  term  becomes: 
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g  dz  dt 
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To  compute  dp/dt,  p  is  substituted  from  the  integrated  form  of  the 

hydrostatic  equation  p  =  p  exp  ( -^z )  where  the  subscript  k  denotes 

*  RT 

a  vertical  level  and  the  overbar  represents  a  vertical  layer  average 
from  k  to  k-1 .  The  partial  derivative  (d/dt)  is  taken  to  obtain  the 
expression: 
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From  Poisson's  equation: 
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where  K  =  — 
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Substituting  Eq.  (6.6)  into  (6.5)  and  decomposing  the  averages 


yields  the  expression  for 
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To  close  this  equation,  the  pressure  tendency  at  the  ground 

(compressible  version)  or  a  time  difference  of  the  perturbation  pressure 

at  the  model  top  (incompressible  version)  is  computed.  These  will  be 

described  below.  The  potential  temperature  tendency  is  computed  from 

dQ.Jdt  of  the  current  time  step, 
ll 

6.1.2  The  Compressible  Hydrostatic  Model 
Aside  from  the  vertical  velocity  diagnosis,  compressibility  can 
enter  a  hydrostatic  model  in  Cartesian  coordinates  through  the  form  of 
the  boundary  condition  for  the  hydrostatic  equation.  This  is  the  case 
for  the  compressible  hydrostatic  version.  The  pressure  boundary 
condition  is  a  prognostic  surface  pressure  derived  from  the  substitution 
of  the  hydrostatic  equation  into  the  compressible  continuity  equation. 

JL  =  9poUi  (6.8) 

3z  dt'  8  dx. 


Integrating  from  the  ground  (z  )  to  the  height  where  p  =  0  (z  )  and 

g  P 

assuming  the  boundary  condition  that  dp/dt  =  0  at  p  =  0. 
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Since  it  is  obviously  not  possible  to  extend  a  Cartesian  model  to  p 


=  0,  a  one-layer  "sigma-p"  model  is  attached  to  the  top  of  the 


HI 


Cartesian  model.  This  provides  the  required  divergence  which  is  assumed 
to  apply  from  the  top  of  the  Cartesian  model  (z^)  to  p  -  0.  The  surface 
pressure  tendency  equation  then  becomes: 
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(6.10) 


where  a  =  p/p#  and  p#  is  the  pressure  at  the  base  of  the  o^  layer. 

The  a  layer  is  configured  as  follows:  The  horizontal  stagger  of 
P 

the  variables  is  identical  to  the  main  model.  In  the  vertical,  p,  and  a 

are  defined  at  a  =1,  while  u,  v,  Q,  and  z  are  defined  at  a  =  .5.  The 

top  boundary  condition  is  a  =  0  at  a  =  0.  The  equations  which  govern 

the  a  layer  are: 
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Note  that  since  p,  is  given  by  the  hydrostatic  integration  from  the 

ground,  the  continuity  equation  for  the  a  model  is  not  required. 

P 

One  of  the  major  advantages  of  the  compressible  model  is  that  mass 
is  implicitly  conserved  since,  in  a  hydrostatic  model,  mass  is  directly 
related  to  pressure.  However,  the  price  that  is  paid  is  in  slower 
execution  time  since  the  time  step  is  limited  by  the  stability  criterion 
for  the  Lamb  wave,  which  travels  at  the  speed  of  sound.  Therefore,  the 
compressible  version  has  been  formulated  in  a  time-split  mode  similar  to 
TC  which  splits  the  Lamb  wave  and  part  of  the  external  gravity  wave  onto 
the  small  time  step.  A  summary  of  the  computational  procedure  for  the 
time  split  mode  follows: 

1)  All  tendency  terms  for  the  prognostic  variables  except  for  the 
horizontal  pressure  and  geopotential  gradients  and  a  30/3o  are  computed 
over  the  interval  t  -  At,  to  t  +  At  . 

2)  The  horizontal  velocities  at  t  -  At,  are  stepped  forward  to 
t  -  At  +  At  with  the  horizontal  pressure  and  geopotential  gradients 

L  S 

and  large  time  step  tendencies.  9  at  o  =  .5  is  stepped  forward  with 
'  39 

0  a  . 
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3)  The  vertical  velocity  is  calculated  and  the  surface  pressure  at 
t  -  At,  +  At  is  evaluated  in  a  backward  step  with  horizontal  velocities 

L  S 

at  t  -  At.  +  At  . 

L  s 


4)  The  pressure  profile  (which  is  constant  over  the  small  time 


m 


steps  at  time  t  since  6  has  no  component  on  the  small  time  step)  is 
shifted  by  the  boundary  condition. 

5)  Steps  2-4  are  repeated  up  to  t  +  At^. 

This  procedure  effectively  isolates  the  Lamb  wave  on  the  small  time 
step  and  also  the  portion  of  the  external  gravity  wave  which  is  produced 
by  the  0^  layer.  The  execution  time  for  a  compressible  hydrostatic  run 
can  be  reduced  by  a  factor  of  2  or  3  by  this  technique. 

6.1.3  The  Incompressible  Hydrostatic  Model 

The  incompressible  hydrostatic  model  uses  the  internal  gravity  wave 
radiation  boundary  condition  of  Klemp  and  Durran  (1983).  This  condition 
is  derived  from  the  linear.  Boussinesq  hydrostatic  irrotational 
equations.  By  transforming  these  equations  to  wave  space  and  allowing 
only  upward  propagating  solutions,  a  boundary  condition  is  derived. 


where  p'  and  w  are  the  Fourier  coefficients  of  p'  and  w,  k  is  the 

horizontal  wavenumber  and  N  is  the  Brunt-Vaisala  frequency.  In 

practice,  the  horizontal  w  field  at  z^  is  Fourier  transformed,  the 

coefficients  are  multiplied  by  Na  /|k|,  and  the  result  is  reverse 

o 

transformed  to  obtain  the  perturbation  pressure  field  at  z^,.  This 
serves  as  the  pressure  boundary  condition  for  the  hydrostatic  equation. 

Although  this  condition  has  been  derived  from  a  simplified  equation 
set,  Klemp  and  Durran  perform  several  scale  analyses  which  show  that 
this  condition  may  still  be  valid  in  the  presence  of  linear, 

rotational,  and  compressible  effects.  Preliminary  numerical  tests  have 
shown  this  to  be  the  case.  Therefore,  this  boundary  condition  can  be 
run  with  any  of  the  three  forms  of  the  vertical  velocity  diagnosis. 
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The  time  step  of  the  incompressible  model  is  limited  by  the  speed 
of  the  external  gravity  wave  which  is  partially  contained  on  the  large 
time  step.  Hence,  there  is  no  benefit  in  the  time-splitting  mode. 
However,  the  incompressible  is  still  nsnally  more  efficient  than  the 
compressible.  The  external  gravity  wave  speed  is  proportional  to  the 
height  of  the  model  top  (z^.) .  In  the  compressible  version  with  the  de¬ 
layer,  the  effective  model  top  is  much  higher  than  z^  since  the 
variables  are  defined  at  a  =  .5.  Even  though  the  terms  of  the  a layer 
responsible  for  the  external  wave  propagation  are  on  the  small  time 
step,  the  external  wave  of  the  main  model,  which  is  partially  on  the 
large  time  step,  manifests  itself  in  the  a layer  throngh  the  dp,/dx 
term.  The  net  result  is  that,  as  z^,  become  larger  (as  the  external  wave 
speed  approaches  that  of  the  Lamb  wave),  the  compressible  version 
efficiency  approaches  that  of  the  incompressible  version. 

6.1.4  Plans  for  Farther  Development 

The  current  plans  for  the  further  development  of  the 
hydrostatic/non-hydrostatic  CSU  cloud/mesoscale  model  include: 

-  initialization  with  observed  data  fields 

-  continued  development  of  a  convective  parameterization 

-  option  for  a  latitude-longitude  grid 

-  inclusion  of  the  long  and  shortwave  radiation  parameterization  of 
Chen  and  Cotton  (1983) 

-  inclusion  of  the  planetary  boundary  layer  and  surface  energy  and 
moisture  budget  parameter izat ions  of  Mahrer  and  Pielke  (1977). 

6 .2  The  Data  Assimilation  and  Analysis  Package 

As  the  first  step  toward  the  initialization  of  the  CSU 
cloud/mesoscale  model  with  observed  data,  a  package  of  data  assimilation 


and  analysis  programs  has  been  written.  These  programs  have  been 
designed  to  be  a  general  and  flexible  data  analysis  system  which  produce 
data  sets  that  can  be  used  for  a  variety  of  purposes  ranging  from 
synoptic  analysis  to  the  initial  and  boundary  conditions  for  numerical 
models  . 

The  package  has  the  ability  to  assimilate  data  from  several 
sources.  The  current  version  accesses  as  the  basic  data  sets  the  NMC 
2.5°  latitude-longitude  mandatory  level  pressure  data,  NMC  mandatory  and 
significant  level  rawinsonde  report,  and  USAF  1°  average  surface 
elevation  data. 

The  assimilation  procedure  is  as  follows. 

-  The  NMC  2.5°  pressure  data  are  accessed.  These  are  originally 
hemispheric  fields  which  are  reduced  down  to  a  grid  that  is  specified  by 
the  user.  This  grid  is  then  stored  as  an  output  file  which  can  be 
plotted  and  analyzed  by  the  plotting  routines. 

-  The  pressure  data  grid  is  interpolated  to  isentropic  surfaces  and 
combined  with  the  surface  elevation  data  to  form  a  complete  synoptic 
scale  2.5°  latitude-longitude  isentropic  data  set  on  the  user-specified 
grid  size.  This  data  set  is  again  stored  as  an  output  file  to  be 
plotted  and  analyzed  if  desired. 

-  The  rawinsonde  reports  are  accessed  and  interpolated  to  the 
isentropic  surfaces.  Options  exist  at  this  point  to  include  only  a 
specified  number  of  the  NMC  reports,  eliminate  any  of  the  NMC  reports, 
or  include  any  special  observations  at  the  analysis  time  that  were  not 
included  in  the  NMC  data  set. 


-  Objective  analysis  are  done  to  the  same  or  different  resolution 
horizontal  latitude-longitude  grid.  This  grid  could  be  a  coarser  grid 
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for  general  circulation  studies,  the  same  2.5°  for  synoptic  studies,  or 
a  finer  resolution  for  mesoscale  analysis  or  model  initialization.  The 
objective  analyses  are  done  on  the  isentropic  surface  using  the  Barnes 
(1973)  technique  utilizing  both  the  rawinsonde  reports  and  the  2.5°  data 
as  observation.  Optionally,  only  the  rawinsonde  reports  can  be  included 
in  the  objective  analyses  if  the  data  coverage  is  adequate. 

-  The  final  product  is  an  "enhanced”  isentropic  data  set  on  a 
latitude-longitude  grid  which  can  be  plotted  and  analyzed  or  used  is  the 
initial  conditions  for  a  numerical  model. 

6 .3  Plans  for  Merger  with  Pielke's  Model 

Recently  an  agreement  has  been  reached  with  Dr.  Roger  Pielke 
wherein  our  modeling  group  and  Dr.  Pielke's  will  jointly  use  the  PP1 
preprocessor  package  and  incorporate  physical  modules  developed  by 
Pielke's  group  as  options  in  the  cloud/mesoscale  model  system.  These 
options  would  include  his  surface  energy  budget  model,  the  ABL 
parameterization,  and  possibly  radiation  schemes. 

This,  we  believe,  will  reduce  our  developmental  time  and  further 
increase  the  flexibility  of  the  model. 


7.0  Investigations  of  Wind  Shears  Produced  by  Precipitating  Convective 


Clouds 

Strong  wind  shears  are  often  associated  with  precipitating 
convective  cloud  downdrafts.  Some  of  the  early  work  in  this  area  was 
based  on  analyses  of  flight  records  (Fujita  and  Byers,  1977;  Fujita  and 
Caracena,  1977)  from  aircraft  which  experienced  sudden  indicated 
airspeed  drops  in  the  vicinity  of  thunderstorms.  More  recent  analyses 
by  Fujita  (1981),  Fujita  and  Wakimoto  (1981),  Mueller  and  Hildebrand 
(1983)  and  Kessinger  et  al .  (1983)  have  disclosed  that  downdrafts  and 
their  associated  near  surface  wind  shears  can  exhibit  substantial 
temporal  variability.  Extreme  wind  shears  were  found  to  occur  in  the 
lowest  few  hundred  meters. 

Similar  variability  in  low  level  shears  have  been  analyzed  in  South 
Park,  Colorado  thunderstorms  with  an  example  shown  in  Fig.  7.1.  In  this 
case  a  very  strong  wind  shear  (~  30  ml s  over  a  500-1000  m  vertical 
distance)  occurred  along  an  inflow/outflow  boundary.  This  shear  zone 
exhibited  significantly  less  magnitude  9  min  later.  Preliminary 
modeling  activities  with  the  3D  cloud/mesoscale  model  have  also  produced 
extreme  transient  wind  shears  at  the  base  of  downdraft  circulations. 

The  example  shown  in  Fig.  7.2  from  model  studies  reveals  significant 
vertical  shear  of  horizontal  winds  in  addition  to  horizontal  variations 
in  horizontal  winds  at  the  base  of  the  downdraft. 

Further  experiments  will  address  the  transience  and  evolution  of 
such  shear  zones  associated  with  downdraft  circulations. 

For  extreme  shear  near  the  surface  (the  lowest  ~  1  km),  adiabatic 
mixed  layers  are  needed.  For  a  given  rainfall  rate,  a  deeper  ABL  will 
in  general  produce  stronger  downdrafts  and  associated  shear.  When  the 
ABL  depth  exceeds  3  km  as  it  often  does  in  dry  climates,  almost  any 
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Figure  7.2:  2D  simuletion  initilized  with  the  26  July  1977  1700  MDT 

South  Park  sounding.  Microphysical  quantities  are  plotted 
in  the  top  panel,  with  the  symbols  representing:  0  - 
rainwater,  x  -  graupel,  +  -  ice  crystals,  and  o  -  cloud 
water.  In  the  lower  panel  perturbation  pressure  contours 
are  super iaq>osed  on  the  relative  flow  field. 
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precipitating  clond  is  capable  of  generating  strong  downdrafts.  For  a 
shallow  ABL,  downdraft  aiagnitndes  tend  to  be  a  function  of  rainfal  rate, 
so  that  only  heavily  precipitating  systems  are  likely  to  produce  strong 
downdrafts  and  wind  shear. 


1 A  1 


8.0  Recent  Changes  in  the  3D  Cloud/Mesoscale  Model  Sent  to  Kirtland 

AFB. 

Since  January  1983,  when  the  CSU  3D  Cloud/Mesoscale  model  was 
installed  at  Kirtland  AFB,  some  major  changes  have  been  made.  These 
change  are  as  follows: 

8.1  Prediction  of  p1 

Because  of  apparent  instabilities  associated  with  prediction  of 
p#  ,  we  made  a  change  to  predict  p'  instead.  In  doing  so,  we  have  used 
the  truncated  form  of  the  pressure  equation  (described  by  Klemp  and 
Wilhelmson,  1978).  Following  the  notation  of  Tripoli  and  Cotton  (1982), 
the  pressure  equation  is  now: 


,  ,  ,  3u.p  0 

op  yp '  t  o  o 

3t  p  0  3x, 

o  o  i 


=  0 


(8.1) 


There  are  no  large  time-step  terms.  The  terms  dropped  were  shown  by 
Klemp  and  Wilhelmson  (1978)  to  have  no  apparent  effect  on  the 
meteorological  dynamics. 

The  w  equation  then  becomes: 
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(8.2) 


The  other  equations  remain  unchanged. 

_  00’  0n ' 

Because  the  part  of  ^  is  evaluated  on  the  large  time  step, 

this  set  of  equations  seems  to  be  capable  of  a  longer  large  time  step, 

0> 

limited  only  by  the  g  —  part  of  the  gravity  wave  acceleration  due  to  w. 

o 

Since  finite  differencing  is  not  involved,  the  gravity  wave  associated 
time  step  limitation  is  from  the  Brunt-Vasal  1  la  frequency  (N)  such  that: 
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At  L  a  ( 


A.  Ma/2 

0  dz1 


(8.3) 


where  a  is  some  constant.  We  think  a  is  between  0.2  and  0.1  since  we 
need  about  5-10  points  in  tine  sampled  as  a  wave  goes  by  to  properly 
resolve  it  without  aliasing. 

8.2  Tod  Boundary  Condition 

We  have  added  the  Klemp  and  Dnrran  (1983)  gravity  wave  radiation 
top  boundary  condition  as  an  option  replacing  the  acoustic  radiation  top 
boundary  condition  option.  This  option  is  implemented  only  for  two- 
dimensional  simulations  thus  far.  The  extension  to  3  dimensions  will  be 
relatively  trivial.  In  essence  the  top  boundary  condition  is 


where  k  is  the  wave  number  and  the  refers  to  the  value  in  wave  space. 
The  condition  was  implemented  on  the  small  time  step  and  hence  a  major 
recoding  of  the  small  time  step  routine  ACODSTC  was  made.  For  the  case 
of  the  boundary  condition  use 

.  AC  X 

.  SE  EX  =  1  . 

In  order  to  implement  this  change  the  NZP  thermodynamic  grid  point  was 
redefined  to  be  coincident  with  the  NZ  w  point.  The  NZP  thermodynamic 
point  is  then  predicted.  In  this  way  a  prediction  of  w*  and  p'  will 
occur  at  the  same  point  so  Eq.  (8.4)  can  be  used.  A  fast  fourier 
transform  routine  is  called  and  will  also  be  needed  in  the  Eirtland 


A 


library. 


8.3  Redefinition  of  DIR2  and  DIR3 


The  definition  of  DIR2  and  DIR3  have  been  redefined  to  be  NXP  and 
NYP  instead  of  NX  and  NY.  NX  and  NY  are  limits  to  which  scalar  points 
are  predicted  in  the  X  and  Y  directions,  respectively.  NXP  and  NYP  are 
NX  +  1  and  NY  +  1  and  DIR2  and  DIR3  are  the  limits  to  which  grid  points 
are  stored.  This  change  means  X  and  Y  I/O  loops  must  now  go  one  point 
longer.  The  change  was  made  to  accomodate  the  introduction  of  a 
specified  variation  in  space  and  time  of  quantities  beyond  the 
integration  grid. 

8 . 4  Introdnction  of  Packing  Input/Output 

The  ability  to  run  a  core  to  packed  core  input/output  routine  on 
previously  core-contained-only-variables  was  added.  The  scheme  allows 
much  large  domain  sizes  to  be  integrated  with  minimal  loss  of  accuracy. 

When  packing  is  activated  (activation  character  "p”).  all 
specified  packing  variables  are  placed  in  memory  packed  by  slabs  in 
ratios  of  one  to  one,  two,  three  or  four.  Each  variable  may  have  a 
different  ratio.  The  variables  are  unpacked  into  a  slab  stencil  as  with 
the  disk/core  I/O  scheme.  Three  or  five  slab  stencils  may  be  specified. 
NCAR  system  integer  packing  routines  are  used  which  attempt  to  limit 
biased  truncation  errors.  The  results  are  impressive.  A  four  to  one 
packing  rate  leads  to  a  13%  increase  in  CPU  usage  with  little  difference 
in  the  predicted  results. 

8.5  Globalizina  of  ADVECT 

In  order  to  simplify  the  advection  routine,  the  flow  for  z,  y  and  z 
advection  was  made  into  a  large  global  containing  the  finite  difference 
scheme  and  the  boundary  condition  logic.  This  shortens  the  advection 
scheme  by  60%  and  makes  the  remaining  logic  more  readable. 
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8.6  Changes  in  MICROPKG 

The  continuing  debugging  and  development  process  has  led  to  many 
changes  in  the  microphysics  package.  The  current  package  is  described 
in  Chapter  4. 
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